



RF Project 761054/711055 
Report No. 384 




BASELINE ESTIMATION FROM SIMULTANEOUS SATELLITE 
LASER TRACKING 

George C. Dedes 

Department of Geodetic Science and Surveying 


(H4SA-CB-182750V BASBLIHB BSTIHATIOH BROH H88-22052 

SIR OLT&H BOOS S&TELLITB L&SEB TBICKIMG (Ohio 
State OniT.i 204 p CSCi. 22 A 

Onclas 

G3/13 0140235 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
Goddard Space Flight Center 
Greenbelt, Maryland 20771 


Grant Mo. MSG 5-265, Supp. No. 10 


October 1987 



The Ohio State University 
Research Foundation 

1314 Kinnear Road 
Columbus, Ohio 43212 


BASELINE ESTIMATION FROM SIMULTANEOUS SATELLITE 


LASER TRACKING 


by 

George C. Dedes 


Report No. 384 


Department of Geodetic Science and Surveying 
The Ohio State University 
Columbus^ Ohio 43210-1247 


October, 1987 



DEDICATION 


To Alexandra 
... who never ceases 
to wonder... 


PREFACE 


This project is under the supervision of Professor Ivan I. Mueller, 
Department of Geodetic Science and Surveying, The Ohio State 
University. The science advisor is Dr. David E. Smith, and the technical 
officer is Dr. Gilbert B. Mead, both at Code 601, Crustal Dynamics 
Project, Space and Earth Sciences Directorate, NASA Goddard Space 
Flight Center, Greenbelt, Maryland 20771. The work is carried out 
under NASA Grant No. NSG 5265, OSU Research Foundation Project No. 
711055. 

This report is an expanded version of a dissertation submitted to 
the Graduate School of The Ohio State University in partisd fulfillment of 
the requirements for the PhD degree. 


iii 



ABSTRACT 


Simultaneous Range Differences (SRD's) to Lageos are obtained by 
dividing the observing stations into pairs with quasi-simultaneous 
observations. For each of those pairs the station with the least number 
of observations is identified, and at its observing epochs interpolated 
ranges for the alternate station are generated. The SRD observables 
are obtained by subtracting the actually observed laser ranges of the 
station having the least number of observations from the interpolated 
ranges of the alternate station. On the basis of these observables 
semidynamic single baseline solutions have been performed. The aim of 
these solutions is to further develop and implement the SRD method in 
the real data environment, to assess its accuracy, its advantages and 
disadvantages as related to the range dynamic mode methods, when the 
baselines are the only parameters of interest. 

Baselines, using simultaneous laser range observations to Lageos, 
have also been estimated through the purely geometric method. These 
baselines formed the standards of comparison in the accuracy 
assessment of the SRD method when compared to that of the range 
dynamic mode methods. On the basis of this comparison it was 

concluded that for baselines of regional extent (i.e., up to 3700 km) the 

SRD method is very effective, efficient and at least as accurate as the 

range dynamic mode methods, and that on the basis of a simple orbital 

modeling and a limited orbit adjustment. 

The SRD method is insensitive to the inconsistencies affecting the 
terrestrial reference frame and simultaneous adjustment of the Earth 
Rotation Parameters (ERP's) is not necessary. Therefore, this method 
offers an inexpensive alternative for projects designed to study regional 
plate tectonic motions. 
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Chapter 1 
INTRODUCTION 

1.1 BASELINE ESTIMATION IN THE DYNAMIC AND SEMIDYNAMIC 
ENVIRONMENT 

In the dynamic environment accurate baseline estimation requires a 
highly sophisticated orbital modeling and a baseline-pass geometry 
leading to near cancellation of the accumulated along-track and 
cross-track orbital errors caused by the erroneous constraints imposed 
on a large number of estimable quantities (Rao, 1973), the recovery of 
which is not possible due to their reduced data sensitivity. In this 
environment proper implementation of the Terrestrial Reference Frame 
(TRF) requires simultaneous recovery of the Earth Rotation Parameters 
(ERP) or utilization of a consistent set of ERPs obtained through a 
separate step. 

Although fulfillment of these requirements makes it possible to 
effectively recover baselines of global and regional extent, it results in 
low temporal resolution of baseline recovery. 

In the semidynamic environment (Section 2.2), and on the basis of 
simultaneous observations, only regional baselines can be recovered with 
an accuracy compatible to that of the observed laser ranges. The 
maximum regional baseline length effectively recovered in this 
environment depends on whether the simultaneous observations collected 
by the baseline end stations are enough to result in a steady state 
response (Section 4.2.1). This, however, is a function of the satellite 
altitude, and for the Lageos satellite the effective regional extent may 
include baselines of up to 3703 km (Section 4.4.2). 

In the semidynamic environment a relatively simplified orbital model 
is required and only the position and the orientation of each of the arcs 
involved is adjusted to "best" fit the available observations (Section 
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2.2.6). Adjustment for the ERP parameters is not necessary since 
proper implementation of the TRP frame is warranted by the use of only 
simultaneous observations. The rekixing data requirements and the 
limited orbit adjustment make it possible to increase the resolution of 
baseline recovery without any loss on the achieved accuracies, and at 
the same time to substantially decrease the required computations, 
thereby making it possible to effectively implement the semidynamic 
methods with limited computer facilities as, for instance, in the personal 
computer (PC) environment. 

The sophistication of the orbital modeling in the semidynamic 
environment can be further simplified by appropriately transforming the 
observed laser ranges to bring them "closer" to the estimable quantities 
being recovered (i.e., baselines). The term "closer" indicates that on 
the basis of the same orbital model the errors eiffecting the computed 
value of the transformed observations, referred to from now on as 
"observables," are smaller than those affecting the computed value of 
the observations themselves. By bringing the observations closer to the 
estimated baselines, the sophistication of the orbital modeling could be 
further reduced if the performed transformation cancels out the errors 
caused by the model simplifications. 

Transformation of the laser range observations to Simultaneous 
Range Difference (SRD) observables brings them closer to the estimated 
baselines (Pavlis, 1982). The potential of using SRD observables to 
estimate baselines was studied with simulated data by Pavlis (1982). The 
results of this study were very promising not only with respect to their 
accuracy but also with respect to their simplicity. 

1.2 SCOPE OF THIS INVESTIGATION 

For the reasons mentioned in the previous section, it was considered 
appropriate and worthwhile to pursue the present investigation, the aim 
of which is to further develop and implement the SRD method (Section 
2.2.1) in the real data environment, to assess its accuracy, its 
advantages and disadvantages as related to range dynamic mode 
methods, when the baselines are the only parameters of interest. 
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Since during the MERIT Main Campaign many stations collected 
simultaneous observations (Section 3.4), it was decided to proceed not 
only with the development of the SRD method but also with baseline 
estimation through the geometric method (Section 4.3). These baselines 
formed the standards of comparison in the accuracy assessment of both 
the SRD and the range dynamic mode methods (Section 4.5). 

In pursuing this study, the geometry does not take part in the 
physical events, as happens in general relativity, but rather it is used 
in its deductive form either purely (i.e., geometric methods) or in 
combination with an inductive form (i.e., dynamic and semidynamic 
methods), both of which were employed to formulate the spatial-temporal 
relationships of the observing stations and the observed satellite 
positions. 

In a pure deductive mode the geometric method is entirely based on 
Euclidean geometry, without any reference to the inductive inference 
that the satellite moves along the path chosen by its physical 
environment. The implied Euclidean geometry is revealed, in the 
arithmetic framework, on the basis of Cartesian coordinates (Section 
2.1.1). There exist, however, configurations in which some of the 
vertices of the resulting figures are free to move along a locus of 
p>oints, thereby forming configurations that are not unique and are 
referred to as "critical configurations" (Section 2.1.3). 

In the inductive mode, the differential form of a satellite’s motion, 
modeled on the basis of Newton’s inductive laws, is expressed through 
its equations of motion (Section 2.2.5). These equations are integrated, 
on the basis of the deductive mathematical framework, to reveal the 
geometric path of the satellite (i.e., orbit), chosen by the satellite’s 
physical environment (Section 2.2.3). Having the geometric path of the 
satellite, the observing stations are connected with this path through 
the Euclidean distance formulated in terms of Cartesian coordinates. 
Since the simplest geometry of the satellite orbit reveals itself in an 
inertial reference frame, the satellite’s equations of motion are expressed 
with respect to such a frame (Section 2.2.4). However, the estimated 
Cartesian coordinates of the observing stations are referenced with 
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respect to an earth-fixed frame (Section 2.2.2). 

The deductive and/or inductive formulation described above contains 
a large number of slow varying quantities which can be considered 
constant for the time span of the observations and a subset of which 
constitutes a set of quantities that eu*e estimable only if the necessary 
units, constants and/or constraints have been properly adopted for 
their unique determination (Rao, 1973). 

There are three types of estimable quantities that can be 
differentiated from their relation to the physical environment, or the 
observing environment, or the links of these two environments. 

In satellite geodesy the interstation distances (i.e., baselines) are 
estimable quantities related to the observing environment. The estimable 
quantities of the physical environment are those related to its cause 
(i.e., Anm and Bnm potential coefficients) and those to its effect (i.e., 
geometric and dynamic characteristics of the orbit). The estimable 
quantities molded in the link of the physical and the observing 
environments are those resulting from the latter as relates to the cause 
and effect duality of the former (i.e., station geocentric distances, 
latitudes, and longitude differences). 

The baselines are computed from the earth-fixed station coordinates 
which are recovered through an inversion process, such as 
least-squares adjustment, on the basis of both the geometric and the 
SRD methods (Sections 2.1.2 and 2.2.6). The input to this inversion 

process are the Simultaneous Ranges (SR) and SRD observables. These 
observables were generated through an interpolation of the observed 
laser ranges because it is quite unlikely, if not impossible, to record 
simultaneous observations to a passive satellite (Section 3.6). Because of 
the peculiarities of the SLR system (Section 3.1 ), it is very likely that 
the recorded laser ranges will be affected not only by white noise but 
also by large blunders (Section 3.2). For this reason and since the SR 
and SRD observables will be generated via an interpolation, it is 
important to edit the laser ranges before proceeding with the geometric 
and the orbit adjustments. 
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In the geometric method the recovered baselines will only be 
affected by the errors resulting from an improper or from the not yet 
reached steady state response (Sections 4.2). The steady state response ' 
of the geometric method, on the basis of a minimum least-squares 
solution, is srffected only by the observational errors. Such a response, 
however, was not possible for longer baselines (Section 4.3.1), and 
therefore an overconstrained solution was adopted to form the standards 
of comparison with the anticipation that it is the least erroneous when it 
is compared either to SRD or to dynamic solutions (Section 4.3.2). 

The accuracy of the baselines, recovered via the SRD (Section 4.4) 
is assessed from the comparison with the baselines obtained via the 
geometric methods and the range dynamic mode methods (Section 4.5). 
The response of the SRD method to the simplification of the orbital 
model has also been investigated (Section 4.6). 
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Chapter 2 

ESTIMATION METHODS 


In this chapter an attempt is made to briefly describe the mathematical 
models and the principles involved in the geometric and dynamic mode 
methods as they are applied to satellite geodesy. Although the 
geometric solutions performed in this study are only used as standards 
of comparison (Sections 4.3 and 4.5), their mathematical formulation is 
presented first because historically the geometric methods were the first 
ones to result in accurate differential positioning. 

2.1 GEOMETRIC METHODS 

Geometric methods are based on the analysis of the relative geometry of 
the observations without any reference to the physical processes 
creating the problem under question. On the contrary, some or all of 
the systematic corrections applied to the observations are computed with 
the use of physical models. 

In the geometric approach of satellite geodesy (Veis, 1960; Mueller, 
1964a), the observed satellite positions are treated as auxiliary 
independent points in space, and they are only used to relate the 
observations geometrically. This in turn leads to the generation of 
space networks. These networks manifest not only the relative geometry 
of the observations, but also any a priori information which is necessary 
for their realization. Thus each observation relates the position of the 
observing station with the observed satellite position. The unknown 
parameter vector includes the Cartesian coordinates of the observing 
stations together with the Cartesian coordinates of the observed satellite 
positions at each of the observing epochs. Since the coordinates of the 
unknown satellite positions constitute an independent set of unknowns, 
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it is necessary to have a sufficient number of observations at each 
observing epoch. The number of observations should be sufficient not 
only to eliminate the unknown satellite coordinates at each of those 

epochs* but also to solve for the unknown station coordinates. 

The process described* however* necessitates the usage of 
simultaneous (referenced to satellite time) observations without any 
reference to the fact that the motion of the satellite obeys the physical 
laws of dynamics. These two distinct features create the advantages 
and disadvantages of the geometric approach to satellite geodesy. 

2.1.1 Mathematical Model 

The mapping of the parameter space into the observational space is 

referred to as observational modeling. The analytical expression 

responsible for the realization of this mapping is referred to as either 
the mathematical or the observational model. The mathematical model 
employed in the geometric approach is that of the Euclidean range from 
a ground station to an observed satellite position expressed in terms of 
Cartesian coordinates: 

Fij = [(Uj-uO^ + (Vj-v<)* + (wj-w<)*]^ - Tij = 0 (2-1) 

The quantity r< j is the true value of the range observable from the 
ground station i to the satellite position j (see Section 3.1)* while the 
quantities Uj, v<* w^ and Uj, Vj* wj denote the true values of the 
Cartesian coordinates of station i and satellite position j. These 
coordinates may be referred to any arbitrarily chosen Cartesian 

coordinate system since ranges are invariant under any rigid body 
rotations. 

The linearized form of equation (2-1) forms the basis for the 
generation of the observation equations when four or more stations are 
observing simultaneously (see Section 2.1.3). With these observation 
equations the normal equations are derived on the basis of a weighted 
constrained least squares adjustment (Uotila* 1967; Krakiwsky et al., 
1967). The resulting normal equations are reduced by eliminating the 
unknown satellite coordinates. The reduced normals are then solved to 
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estimate the stations' Cartesian coordinates which are finally transformed 
to interstation distances. 

Inverting the normals on the basis of the minimum required 
information (i.e., minimum constraints) leads to baseline errors that 
depend solely on the errors of the observed ranges and on their 
geometric strength as well. This is true only when the scale has been 
properly incorporated into the solution either through the observations 
and their geometry, or if this is not possible, through some other 
additional constraints (see Section 4.3). 


2.1.2 Normal Equations 

The observation equations used to derive the normal equations are 
obtained from the linearization of equation (2-1). Linearization is 
achieved with a Taylor series expansion about the approximate values of 
the station coordinates, the satellite coordinates and the observed 
ranges as well. The resulting linear equations have the form: 


^ij “ V,j + Lij - 0 

where 

Aij = [ aij : “^ij] 


3 (Uj, Vj, Wj) 


( 2 - 2 ) 


(2-3) 

(2-4) 



du< 
dvj 
dw. 
du{ 
dv< 
dw, j 


(2-5) 


V<j = residual vector corresponding to the range observable 
L< j = the computed minus the observed range 
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In cases when the geometric strength of the observations is not 
good enough to warrant a steady state response, it could still be 
possible to reach such a response if estimated or observed interstation 
distances are incorporated into the solution (see Section 4.3). This is 
accomplished by introducing the inter station distances as fictitious 
observations into the adjustment. Appropriate weights should be 
applied to these (fictitious) observations in order to avoid any scale 
conflicts that might contaminate the solution (Uotila, 1967). 

The mathematical model of the interstation distance between stations 
k and i has the following form 

Gfci = [(ui-Uk)* + (vi-Vk)2 + (wi-Wk)^]’^ - Lki = 0 (2-6) 

where Lki is the true value of the fictitiously observed range between 
these two stations. Linearization of equation (2-6) about the 
approximate station coordinates and the fictitiously observed distance 
results in 

Cki Xki - Vki + Dki = 0 
where 

Cki = [ Ti ; -Ti ] 

3Gki 

T A = -- - - 

a(ui, Vi, wi) 

Vki = residual of the (fictitiously) observed interstation distance 

Dki = the computed minus (fictitiously) observed interstation 
distance 

So far the observation equations have been developed on the basis 
of one range observation and one interstation distance. Considering 
many range observations and many fictitiously observed interstation 
distances, equations (2-2) and (2-7) take the following form: 


(2-7) 

( 2 - 8 ) 

(2-9) 
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( 2 - 10 ) 

( 2 - 11 ) 


AX - V + L = 0 
CX - Vc + D = 0 

These equations can be rewritten as 


■ A ■ 


■ L ■ 


■ V 


X + 


— 


C 


_ D _ 


. Vc . 


( 2 - 12 ) 


and with some obvious substitutions they take the following form: 

A*X + L* = V* (2-13) 

Equation (2-13) forms a set of observation equations which are used to 
derive the normal equations on either deterministic or statistical 
grounds. 

Deterministically the principle of least squares requires that the 
quantity (v*Tp*v* + X^P^X) assumes minimum value, subject to the 
condition A*X - V* + L* = 0. The matrix Px is the weight matrix 
associated with the coordinates of the ground stations and of the 
observed satellite positions, while the weight matrix P* takes the 
following form 



(2-14) 


where P and Pc are the weight matrices associated with the range 
observables and the (fictitiously) observed interstation distances 
respectively. 

Statistically the maximum likelihood principle requires maximization of 
the a posteriori conditional density function of the parameter vector X 
given that the observations L* have been obtained. The resulting 

A 

estimator X is referred to as Bayesian least squares estimator. Both of 
the above principles lead to the same estimator X only if normality is 
assumed not only for the a priori density function of the estimated 
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parameter vector X, but also for the conditional density function of the 
observational vector L* given that the parameter vector X has been 
estimated. These assumptions should only hold for the maximum 
likelihood principle since any least squares estimator is a distribution 
free estimator. The estimator X is obtained from the resulting normal 
equations (Uotila, 1967; Krakiwsky et al., 1967; Cappellari et al., 1976): 


ATPA + Px 
C 


CT 



= 0 


(2-15) 


where kc are the Lagrange multipliers associated with the (fictitiously) 
observed interstation distances. Elimination of the Lagrange multipliers 
from equation (2-15) leads to the following form for the estimator X: 

X = -(ATPA + CTPcC + Px)-‘ (ATPL + CTPcD) (2-16) 

Substitution of equations (2-2) through (2-5) and (2-7) through 

(2-9) into equation (2-16) followed by elimination of all the satellite 

coordinates leads to the following equations (Krakiwsky et al., 1967): 

- The 3x3 diagonal matrix associated with the k^^' ground station 

Nkk = (Z axjT P^j a^j) - I {akjT p^j a^j [z a^j^ p.^ a<j] 

Pkj akj} + Pk + [ J ] TfcT p^^ Tk (2-17) 

- The 3x3 off-diagonal matrix corresponding to the k^^ and ground 
station 

Nki = -I (akj^ Pkj axj [z a,jT p^j a,j] a^jT PijaijJ + 

[ J ] TfeT Ti (2-18) 

- The 3x1 constant vector associated with the k^^ ground station 
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(2-19) 


U|c - - (Z Pkj Lfej] + Z (akj^ Pkj Skj [| Pij aij) 

J J 

Z a^jT P^J L,j } + [ 5 ] TJ p^t 

In equations (2-17) and (2-19) the j summation is performed over all 
the satellite positions observed from station k, while in equation (2-18) 
the j summation is performed over all the satellite positions observed 
simultaneously from both stations k and i respectively. The summation 
i, on the other hand, is performed over all the stations observing the 
satellite position j simultaneously. In these three equations the number 
1 is used only when the interstation distance between the stations k 
and i is involved in the solution, otherwise the value of 0 is used. 

Equations (2-17) through (2-19) form the basis for the estimation of 
the ground station coordinates by inverting the normals through a 
procedure accredited to Banachiewich. This procedure is carried out in 
two steps. 

- The first step involves the representation of the normal matrix as 
a product of right and left triangular matrices with the left 
triangular matrix having unit diagonal elements. 

- The second step involves the computation of the inverse normal 
matrix on the basis of only the inverted diagonal elements of the 
right triangular matrix and the off-diagonal elements of the left 
triangular matrix. 

The above procedure is very closely related to the Cholesky algorithm 
(Uotila, 1967). 

2.1.3 Critical Configurations 

In the geometric approach of satellite geodesy the observed satellite 
positions (targets) are treated as auxiliary independent points in space. 
They are used only to relate the observations geometrically through the 
resulting range space networks (see Section 2.1). In certain cases the 
ground stations and/or the targets which form the vertices of this 
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network are involved in a kind of configuration for which a unique 
adjustment is impossible although the number of observations is 
sufficient and the coordinate system well defined. These configurations 
au*e referred to as "critical." 

With range observations the critical configurations have been 
extensively studied in the past (Rinner, 1966; Blaha, 1971; Tsimis, 1972, 
1973). The critical configurations have been traditionally analyzed 
according to whether all of the observing stations are either in a plane 
or generally distributed in space. For both of these cases the resulting 
singularities are divided into three categories: 

1. Singularity A, resulting from the relative geometry of an individual 
station connected to its observed targets. 

2. Singularity B, resulting from the relative geometry of the observing 
stations only. 

3. Singularity C, resulting from the relative geometry of all the 
observing stations connected to their observed targets when 
singularity A and singularity B are not present. 


2.1.3.1 Critical configurations when all of the observing stations 
are in a plane . When all of the observing stations are in a plane the 
singularity problem is analyzed according to the number of stations 
observing all the targets. This number may be three or more, or less 
than three. The number three is important since ranges from three 
stations are needed to eliminate the coordinates of one target, provided 
that this target is not located on the plane of these three stations 
(Blaha, 1971). 

If the number of ground stations observing all the targets is three 
or more, singularity A occurs when an individual station — excluding the 
stations used to define the coordinate system, since for these three 
stations singularity A cannot occur— is either observing less than three 
distinct targets or is in the same plane with all of its observed targets. 
Furthermore, singularity B occurs when all the observing stations but 
one are lying in a straight line or more generally when all the stations 
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are lying on a second-order curve. Since at least five stations are 
needed to determine a second-order curve, singularity B can only be 
avoided if six or more stations are involved. In the absence of 
singularities A and B singularity C occurs when the stations making 
off-plane observations (i.e., observed targets are not in the same plane) 
are not themselves off-curve stations (i.e., not lying in the same 
second-order curve). To avoid singularity C at least three off-curve 
stations should make off-plane observations (Blaha, 1971; Mueller et al., 
1975). In the case when all the stations observe all the targets, 

singularity A loses its importance because it always implies singularity 
C, since off-plane observations are necessary to avoid singularity C. 

When there do not exist three stations observing all the targets, 
elimination of the coordinates for all the targets using the same three 
stations cannot be achieved. Thus, in the elimination process one, two 
or all of these three stations will have to be replaced. This leads to 
the first, second and third replacement respectively, and therefore to at 
least four-station events. 

We first denote with k the station used in the first replacement. In 

this replacement singularity A for all the stations but k, or singularity 

B for all stations would occur as though there were three stations 

observing all the targets. For station k, however, singularity A occurs 
if any new stations coming into play are lying either on the x-axis of 
the local coordinate system (e.g., line formed by two of the three 
stations used to define the local coordinate system) or in the 

intersection line (denoted by i) of the plane tt (see below) with the 
plane of the ground stations. The plane n, if it exists, is the plane 
containing station k together with the satellite positions (denoted by j^) 
that were observed by this station (e.g., k) up to the epoch at which 
the new station(s) started observing. Singularity C, in the absence of 
singularities A and B, is further analyzed according to whether are 
off-plane or in-plane targets. If are off-plane targets, singularity C 
occurs whenever stations making off-plane observations are not 
themselves off-curve stations. The case when j|< form in-plane targets 
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is not discussed here because it is very unlikely to encounter in 
practice. This case, of rather academic interest, is discussed in (Blaha, 
1971, page 63). 

Next we denote by s' the station used in the second replacement. 
In this replacement, if there is no other station s" to start observing 
for the first time after s' has started, singularity A occurs as though 
only the first replacement would have taken place. However, with 
station s" present, singularity A for station k occurs if, in addition to 
the above, the station s" is lying either in the x-axis of the local 
coordinate system or in the line i. Furthermore, singularity A for 
station s' occurs if in addition to being in the plane n' (defined below), 
the station s " is also lying either in the line defined by the station 
used as the origin of the local coordinate system and the station k or in 
the intersection line (denoted by i') of the plane n' with the plane of 
the ground stations. The plane w', if it exists, is the plane containing 
the station s' together with the targets (denoted by jg') observed by 
the station s' up to the epoch at which the station s" started 
observing. Singularity B occurs as though three stations observing all 
the targets exist. If jg' and are off-plane targets, in the absence of 
singularity A and B, singularity C occurs, when no other stations 
besides k and s' exist to make off-plane observations. The case when 
Js' or jj, form in-plane targets is not discussed here because of the 
unlikelihood of encountering it in practice. 

The singularities resulting from three or more replacements are 
similar to the ones described above. By avoiding singularities A, B and 
C a nonsingular network can only be formed if at least six stations in at 
least four station events are involved. This is so because five stations 
are needed to define a second-order curve and only the sixth station is 
possible to serve as an off-curve station. 

Once a nonsingular network has been realized any extension of it 
will result in a nonsingular one if for any additional station singularity 
A is eliminated and if no target is lying in the plane of the ground 
stations. 
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2. 1.3.2 Critical configurations when the observing stations are 
generally distributed in space . When the ground stations are in general 
configuration, singularity B loses its meaning because the effect of 
ground stations cannot be separated from that of the targets. 
Consequently, singularity B will not be considered here. However, 
another type of singularity called "reverse singularity B" is the 
singularity B if one assumes that the satellite points (targets) observe 
the ground stations. Therefore, this singularity occurs when all the 
targets are in a plane in a second-order curve. This in practice could 
approximately occur when two short passes of about the same altitude 
have been observed. 

Having the observing stations in a general configuration, a 
nonsingular network can be formed if at least six targets are being 
co-observed by at least four stations. Accordingly, the analysis of 
critical configurations proceeds by grouping the ground stations in 
tetrads ("quads"). With four stations observing all the targets, 
singularity A occurs only with respect to the fourth station because 
singularity A never occurs for the three stations that have been used 
to define the local coordinate system (Blaha, 1971). With respect to the 
fourth station, singularity A occurs if edl the targets are lying on a 
plane through this station, or if all the targets are on the plane formed 
by the stations used to define the local coordinate system. 
Furthermore, in the absence of singularity A, singularity C occurs when 
all the observing stations and all the targets are lying on a 
second-order surface. 

With more than four stations observing, the singularity problem is 
analyzed by grouping the observing stations in quads. If the number 
of stations observing all the targets is three or more, singularity A 
occurs, as though all the observing stations were lying on a plane (see 
Section 2.1.3.1). In the absence of singularity A, singularity C occurs 
either when all the observing stations with their observed targets are 
located on a second-order surface or when each tetrad of stations 
together with its observed targets are located on specific second-order 
critical surfaces. These surfaces intersect each other in one 
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second-order curve containing the three stations used to eliminate the 
coordinates of each target and to define the local coordinate system. 
Furthermore, when edl the stations are co-observing, these second-order 
critical surfaces coincide. 

When three stations observing all the targets do not exist, then the 
concept of station replacement should be utilized in exactly the same 
way as described in the previous section. Proceeding with this concept 
it is found that singularity A occurs as though three stations observing 
all the targets exist. As for singularity C, it is again associated with 
other specific second-order surfaces in addition to the ones resulting in 
singularity C when three stations observing adl the targets exist. 

By avoiding singularities A and C and reverse singularity B a 
nonsingular network can be formed. What is important to keep in mind 
is that when the ground stations are generally distributed in space a 
nonsingulEU* network can be formed if at least six targets are 
co-observed by at least four stations . In fact four stations and five 

targets can uniquely define a second-order surface, and the sixth target 
could make the network nonsingular if it is not located on this surface. 

Once a nonsingular network has been realized, an extension of it 
will result in a nonsingular one if singularity A is eliminated for any 
additional station and if no targets are on the plane of the three 
stations defining the Cartesian coordinate system. 


2.2 DYNAMIC AND SEMIDYNAMIC MODE METHODS 

In contrast to the geometric methods, the observed satellite positions in 
the dynamic and in the semidynamic methods are not treated as auxiliary 
independent points in space but rather they are constrained to lie in a 
space curve (Schwarz, 1969). This curve should resemble within the 
required degree of accuracy the satellite orbit under question. The 
satellite orbits, on the other hand, are modeled either empirically or 
dynamically or by combining both empirical and dynamical modeling. 

Empirical modeling of satellite orbits was extensively used in the 
early years of satellite geodesy since many of the model parameters 
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entering the equations of motion were not known with the degree of 
accuracy required for precise geodetic work (Mueller, 1964b). These 
methods are used today in very special circumstances and only in 
combination with dynamic modeling (Tapley et al., 1985a). 

Dynamic modeling results in three second-order differential 
equations or six first-order differential equations referred to as 
equations of motion of the satellite. These differential equations are 
integrated either analytically (i.e., general perturbation methods) or 
numerically (i.e., special perturbation methods) to generate the satellite 
orbits. 

In the general perturbation methods, the equations of motion of the 
satellite are reformulated in terms of a set of orbital elements leading to 
a set of differential equations which can be integrated analytically. 
Unfortunately, a closed form analytical solution for the equations of 
motion of the perturbed two-body problem does not exist. It is 
possible, however, to obtain approximate analytical solutions either by 
restricting the complexity of the perturbation model or by truncating 
high power expansions (Kaula, 1961, 1966; Mueller, 1964b; Goad, 1977; 

etc.) These solutions are approximate and in many cases cannot be used 
for precise geodetic work. They are extremely useful, however, in order 
to gain a keener insight into the effects of various perturbing forces on 
the satellite orbits. 

In the special perturbation methods, the equations of motion of the 
satellite are integrated numerically (see Section 2.2.3). The main 

advantage of these methods is that all the perturbing forces can be 
accommodated to a high degree of accuracy. The special perturbation 
methods, on the other hand, have proven to be computationally more 
efficient, if one takes into consideration the high repetition rate of 
recent geodetic observations (Rizos et ed., 1985; Krakiwsky et al., 1985). 

A combination of empirical and dynamic modeling is usually employed 
either when the satellite orbits are integrated continuously over long 

periods of time (i.e., two months or more) or when unexplained 

perturbing accelerations are present. In the latter case the dynamic 
models are supplemented with empirical models for the as yet not fully 
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understood perturbing accelerations while in the former the empirical 
models are employed to account for the accumulated along-track, 
cross-track and radial errors (Tapley et al., 1985a). For instance, the 
draglike acceleration of Lageos’ orbit which causes a decay of the 
semimajor axis at a rate of 1 mm/d is modeled empirically (see Section 
2.2.5). 

The dynamic mode methods are further subdivided into semidynamic 
(short-arc) and dynamic (long-arc) methods. There is no clear 
distinction between these two terms and their exact meaning depends on 
the investigator and on the kind of problem being analyzed. 

In the present study the estimation of the baselines is performed in 
the semidynamic mode environment. In this environment the lengths of 
the arcs employed are relatively short (i.e., mostly up to three days and 
very rarely up to seven days) (see Section 4.4). A relatively short arc 
is defined as having a maximum length over which the total modeling 
error of a simple dynamic model is well below the noise level of the 
observations (i.e., an order of magnitude or less). Consequently, with 
this definition one may select a relatively simple dynamic model and then 
determine the length of the arc, or one may choose the length of the 
arc and then determine the required sophistication for the dynamic 
models. With such a procedure the systematic errors caused by model 
imperfections cannot accumulate up to a level that may corrupt the 
semidynamic mode solution. The relatively short arcs, however, are not 
stable in the sense that their position and orientation in space depends 
primarily on the geometry of the observations. This instability may also 
cause ill-conditioning of the normal equations (see Section 4.4). 
Furthermore, relatively short arcs cannot be well tracked to bring 
tracking sites of a global extent into a consistent satellite reference 
frame. This implies that it is not possible to use semidynamic mode 
methods for absolute position determination. Instead these methods can 
be effectively used for baseline determinations (Latimer et al., 1977; 
Christodoulidis et al., 1981; Pavlis, 1982; Section 4.4). Baseline estimates 
are even more accurate if the observables are insensitive to the 
position, orientation, and the shape of the trajectory as is the case, for 
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instance, with the Simultaneous Range Difference observables (see 
Section 4.4; Pavlis, 1982). With range observables, which are sensitive 
not only to the position and orientation of the trajectory in space but 
also to its shape, it is still possible to obtain accurate baselines if a 
local support tracking network is available (Christodoulidis et al., 1981). 

2.2.1 Simultaneous Range-Difference Semidynamic Mode Method 

On the basis of the discussion presented in the previous section and 
keeping in mind that our aim is to achieve highly accurate differential 
positioning, wc have chosen to use in this investigation the semidynamic 
(short-arc) method formulated in the context of special perturbations 
(see Section 2.2.3). Furthermore, the laser range observations have 
been transformed to Simultaneous Range-Differences (SRDs), and 
although differencing is a noise generating operation it is anticipated 
that these observables are less affected by the biases in the orbit, the 
reference frame and the observations themselves (Pavlis, 1982). 

Using laser range observations to Lageos it is impossible to obtain 
strictly simultaneous observations not only because Lageos is a passive 
satellite but also because there will always exist synchronization errors 
among the various observing stations. Therefore Simultaneous Range 
Differences can only be obtained through an interpolation (see Chapter 
3). More specifically, the observing stations are divided into pairs of 
simultaneously observing stations. For each pair the station with the 
most observations is interpolated to generate ranges at the observed 
epochs of the alternate station. Finally the interstation distances for 
each of the pairs involved are estimated by processing the generated 
SRDs through a least squares adjustment formulated in the context of 
the special perturbation methods as applied in the semidynamic mode 
environment. 

2.2.2 Mathematical Modeling 

The mathematical model for the Simultaneous Range Difference (SRD) 
observable flpj is obtained by subtracting the Euclidean ranges from 
station 2 and station 1 to the simultaneously observed satellite position j 
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SfS = {(sj - ii.)’(sj - x.]f - {(sj - x,)'(sj - x.))* 


( 2 - 20 ) 


where the vectors Sj = (uj, Vj, wj)T, Xi = (ui, Vi, w,)T and Xj = (ua, 
Va» Wa)^ denote the Cartesian coordinates of the satellite position j, 
ground station 1 and ground station 2 respectively. Since the SRDs are 
invariant under any rigid body rotation the above vectors could be 
expressed in any arbitrarily chosen Cartesian coordinate system. In the 
present study, the vectors Sj, Xi and Xa at epoch j are expressed in a 
Cartesian coordinate system whose origin is conveniently chosen to 
coincide with the center of mass of the earth, and its orientation is 
aligned to that of the true-of-date system (Mueller, 1969). 

The adjusted parameters, in any estimation procedure involving a 
dynamic process, are transformed to a reference frame in which they 
can be considered constant for a certain period of time. This period 
should be long enough to allow for collection of a sufficient number of 
observations needed for a reliable recovery of the adjusted parameters. 
For this reason the ground station coordinates are transformed to a 

terrestrial reference frame (TRF) while the coordinates of the satellite at 
epoch j 8u*e transformed to a celestial reference frame (CRF) with the 
help of the following formulas {see Section 2.2.4) 

Sj = NPRj ; X, = sty, , i = 1, 2 (2-21) 

The quantities S, N and P designate the earth rotation, the nutation and 
precession matrices respectively, while the vectors Rj, Yi and Y 2 denote 
the inertial position vector of the satellite at the epoch j and the 

earth-fixed position vectors of stations 1 and 2 respectively. The 

inertial and earth-fixed frames correspond to the CRF and TRF frames 
respectively as they are described in Section (2.2.4). Substituting 

equation (2-21) into (2-20) one obtains 



where 
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Dj^ = NPHj - S^Ya and Dj^ = NPRj - STYi (2-23) 

The satellite position vector Rj is a function of an initial state 
vector and a large number of parameters affecting the motion of the 
satellite (i.e., potential coefficients, reflectivity, etc.)* The choice of the 
model parameters to be estimated depends on the data coverage and 
distribution which in turn dictates the adopted lengths for continuous 
integration of the satellite orbit (see Section 2.2.1). In the present 
study, the shape of each of the satellite arcs involved is assumed known 
and only its position and orientation in space is adjusted to "best" fit 
the available data (see Section 4.4). Thus, the only adjusted parameters 
inherent to the satellite position vector Rj are the components of the 
initial state vector of the corresponding arc. 

In the derivation of the observation equations, on the basis of the 
equation (2-20), one needs the satellite position vectors at each of the 
observing epochs together with their partial derivatives with respect to 
the corresponding initial state vector. The former is obtained by 
integrating the equations of motion of the satellite while the latter is 
obtained by integrating the variational equations of state (see Section 
2.2.3). The partial derivatives with respect to earth-fixed station 
coordinates, also needed in the derivation of the observation equations, 
are easily obtained by differentiating equation (2-20) (see also Pavlis, 
1982). The resulting observation equations are used to obtain the 

normal equations through a weighted least-squares adjustment (see 
Section 2.2.6). The normal equations are subsequently solved to estimate 
the initial state vectors for each of the eircs involved together with the 
earth-fixed coordinates which are finally transformed to interstation 
distances. 

In the present study the initial state vectors are treated as 
"nuisance" parameters, and therefore one is not concerned with how well 
each of those initial state vectors has been recovered as long as the a 
posteriori variance of unit weight is close to unity. In fact, the reason 
for using SRDs instead of ranges is to reduce the need for accurate 
knowledge of the satellite orbits and yet to increase the potential for 
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baseline estimation with an accuracy compatible to or even better than 
that of the observations. This is possible because SRDs have the 
potential to reduce the effects of biases caused not only by the orbital 
model and the reference frames but also by eliminating uncorrectable 
systematic errors affecting the laser range observations (Pavlis, 1982). 


2.2.3 Orbit Determination with the Method of Special Perturbations 
Dynamic and semidynamic methodsi based on special perturbations, 
require integration of the satellite’s equations of motion. The degree of 
sophistication in formulating these equations depends on the integration 
length and the required accuracy. 

Following the MERIT standards the relativistic perturbations are 
ignored from the equations of motion (Melbourne et al., 1983). 
Accordingly, ephemeris time (t) constitutes the independent variable in 
the equations of motion. Up to 1983, ephemeris time was used as an 
independent variable in the planetery equations of motion and therefore 
in the construction of all the almanacs. Since January O’* 1984, 

ephemeris time has been replaced by Terrestrial Dynamical Time (TDT) 
and Barycentric Dynamical Time (TDB) (The Astronomical Almanac, 1984). 
This was a necessity since data collected in interplanetary missions are 
routinely processed in the context of the relativity theory (Moyer, 1971). 
In this context TDT time corresponds to proper time (i.e., time measured 
by the observer’s clock) while TDB corresponds to coordinate time (i.e., 
time measured at the barycenter of a motionless solar system in the 
absence of all gravitational fields). 

At each observing epoch j, the ephemeris time (tj) is computed from 
the Universal Coordinated Time (UTCj) with the help of the following 
formula 

tj = 32?184 + [taI - UTC]j + UTCj = TDT (after 1984) (2-24) 

where 

tj = ephemeris time at the epoch j 
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[TAI - UTC] j = no. of leap seconds at the epoch j (The Astronomical 
Almanac, p. B5) 

UTCj = Universal Coordinated Time at the epoch j. 

Using the ephemeris time (t) as an independent variable, the Lageos 
equations of motion take the following form (Cappellari et al., 1976; 
Pavlis, 1982): 

R - RpM RnS ^TD ®SR ®AT (2—25) 

Each of these accelerations is expressed relative to the center of mass 
of the earth. More specifically 

- d* f ) T 

R = I^Uj, Vj, WjJ = total quasi-inertial Lageos acceleration at 

the epoch j 

RpM = gravitational acceleration due to point masses 

Rms = gravitational acceleration due to nonsphericity of the gravi- 
tational potential 

Rxo = acceleration due to solid earth tidal effects 
Rsr = acceleration due to solar radiation pressure 
R/^t = Lageos empirical drag-like acceleration 

The acceleration vector R is a function of an initial state vector and a 
large number of parameters affecting the motion of the satellite. These 
parameters pertain to the gravitational potential, to solar radiation 
pressure, etc. As it was described in the previous section, the only 
adjusted orbital parameters considered in this study sure the initial state 
vectors of all the arcs involved. With these orbital parameters the 
variational equations assume a very simple form (Pavlis, 1982): 

Y(t) = A(t) • Y(t) (2-26) 
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with initial conditions 

Yo « [ I I 0 ]3x6 (2-27) 

where 


Y(t) = 


3R(t) 

a(R(to), R(to)) 


3X6 


and 


(2-28) 


A(t) 


dR 


3R(t) 


3X3 


(2-29) 


The matrix Y(t) is referred to as the state transition matrix and is used 
to map the variations of the initial state into variations of the current 
state. 

Equations (2-25) and (2-26) can be integrated either by one-step or 
by multi-step methods. In each integration step, the multistep methods 
require fewer derivative evaluations than the one-step methods of 
compatible accuracy. Fewer computations, on the other hand, not only 
reduce the round-off errors but also require less computing time. 
Furthermore, since these methods possess a larger number of parasitic 
solutions they are more susceptible to instability problems. 

The multistep algorithm used in our study employs a self starting, 
variable step, variable order predictor-corrector mode of operation. 
This mode selects the order automatically while the stepsize is subject 
to accuracy requirements and numerical stability. Keeping the stepsize 
constant, the predictor-corrector is reduced to an Adams-Bashforth 
predictor of order q and to an Adams-Moulton corrector of order q + 1. 
With this algorithm the second-order differential equations are 
integrated directly without reducing them to a first-order system, 
because a second-order set exhibits better numerical stability 
characteristics. The described algorithm was developed and implemented 
in computer coded form by Krogh (1969a, 1969b, 1973a, 1973b, 1974). 
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2.2.4 Reference Frames and Systems 

Reference frames constitute realizations of reference systems. The 
reference frames are used to describe the spatial relationships and the 
temporal variations of objects on the Earth (i.e., terrestrial frames) and 
in space including the Earth (i.e., celestial frames) (Kovalesky and 
Mueller, 1981). A reference system consists of an underlying principle 
and all those elements (e.g>, physical environment, theories and 
constants) that are necessary to accomplish its realization. The elements 
of a system, depending on the application and the accuracy 
requirements, are selectively chosen and therefore the term 
"conventional system" is often used to designate the selection process 
that is usually involved in any realization of a reference system. 

In this context the lAU/IUGG MERIT and COTES Joint Working 
Group recommended the following concepts in regard to reference 
systems and frames (Wilkins and Mueller, 1986): 

The Conventional Terrestrial Reference System (CTRS) be 
defined by a set of designated reference stations, theories and 
constants [necessary elements], chosen so that there is no net 
rotation or translation between the reference frame and the 
surface of the earth [underlying principle]. The frame is to be 
realized by a set of positions and motions of the designated 
reference stations. 

The Conventional Celestial Reference System (CCRS) be defined 
by a set of designated extragalactic radio sources, theories and 
constants [necessary elements], chosen so there is no net rotation 
between the reference frame and the set of the radio sources 
[underlying principle]. The frame is to be defined by the 
positions and motions of the designated radio sources. The origin 
of the frame is to be the barycentre of the solar system. 


The above concepts are to be incorporated in the operation of the new 
International Earth Rotation Service (lERS). This service is scheduled 
to start operating as of January 0^*, 1988 (Mueller and Wilkins, 1986). 

In the Newtonian framework, the reference frame implied by a CCRS 
can be considered as being an ideal inertial frame in the sense that the 
time is homogeneous and the space described by this frame is 
homogeneous and isotropic (Landau et al., 1960). In the general 
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relativistic framework, on the other hand, the reference frame realizing 
a OCRS is aimed to describe the curved space time for which a global 
inertial reference frame does not exist (Moritz, 1979; Fukushima, 1986)! 
This is the reason why, in the above recommendations, the term 
"inertial" has been dropped and the term "celestial" has been used 
instead. 

For precise geodetic work, these seemingly conceptual differences 
manifest themselves when 10“* or 10“* accuracies are sought. 
Therefore, when working at such accuracy levels, care should be taken 
to account for relativistic effects either by using Newtonian formalism 
with small corrections (Moritz, 1979) or by formulating the problem 
entirely in the general relativistic framework (Fukushima, 1986). In the 
present study, in accordance with the MERIT standards and since the 
obtained accuracies hardly reach the 10“® level, we have used the 
Newtonian formalism to formulate the equations of motion of the satellite. 

The status today in terms of reference systems and frames is 
confusing because the user community employs a variety of different 
celestial systems (i.e., extragalactic radio source systems, stellar 
systems, dynamical systems, etc.) and a variety of different terrestrial 
systems as well (i.e., BIH terrestrial reference system, CSR terrestrial 
reference systems, etc.). Investigations, however, are currently 
underway with the objective of linking all of the available terrestrial 
systems into a unified terrestrial system referred to as "BIH Terrestrial 
Reference System" (BTS) (Boucher and Altamimi, 1985, 1986). Linking 
different celestial systems, through their frames, into an ideal celestial 
frame is not an easy task not only because of lack of collocations but 
also because daily polar motion resolution is necessary (Mueller, 1985). 
This kind of resolution is not achievable by the satellite related systems 
due mainly to the deficiencies in nutation theory (Himwich and Harder, 
1986) and to inadequate observational coverage. 

The choice of ideal terrestrial and celestial frames is not important 
in baseline estimation. It is important, however, to consistently link the 
involved terrestrial and celestial reference frames, by choosing the 
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appropriate set of transformation parameters. Effective choice of the 
transformation parameters would only assure a reliable recovery of the 
relative geometry of the observations since these parameters cannot be 
effectively recovered in a semidynamic mode environment. The relative 
geometry of the observations manifests the way the satellite arcs are 
related to the observing stations. Reliable recovery of the relative 
geometry, on the other hand, results in accurate baseline estimation 
simply because baselines are estimable quantities. Estimable quantities 
are molded by the geometric and dynamic characteristics creating the 
problem under question. 

In the present study, the Terrestrial Reference Frame (TRF) is 
implied by the gravity field used to integrate the equations of motion 
and by the adopted polar motion series. The origin of the TRF frame, 
relative to the center of mass of the earth, is defined by the potential 
coefficients Cto* Cn and Sn, while its orientation is primarily defined 
by the potential coefficients Cj, and Sai as well as C22 and S22. More 
specifically, the coefficients C2t and S21 define the orientation of the 
third axis, while 022 and S22 define the orientation of the first axis. 
The orientation of the first axis, however, is weakly defined because the 
Earth’s equatorial moments are nearly equal (i.e., C22 » S22)« 

The modified GEML 2 gravity field, proposed by the MERIT standards, 
has been replaced in this study by the PGS 1680 gravity field 
(Christodoulidis et al., 1985 ). This violation of the MERIT standards was 
necessary to make the gravity field consistent with the adopted BIH 
ix>lar motion series. The modified GEML 2 gravity field has all its 
coefficients but S21 and C21 equal to the coefficients of the GEML 2 
gravity field (Lerch et al., 1985 ). C^t and S21 have been modified to be 
consistent with the mean BIH polar motion values, computed over a 
complete wobble cycle which lasts from 6.5 to 7 years (Melbourne et al., 
1983 ). Modifying only the C21 and S21 has caused inconsistencies as to 
what frame the coefficients of the modified GEML 2 field refer. As a 
result, the PGS 1680 gravity field was developed in order to avoid these 
inconsistencies and the resulting confusion as well. In the development 
of this field the coefficients 821 and C21 were constrained to the BIH 
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implied values while the remaining coefficients were free adjusting 
(Christodoulidis et al., 1985). 

The coefficients Ciot Cn and Sn of the PGS1680 gravity field are 
zero thereby imposing the origin of the TRF frame to coincide with the 
center of mass of the earth. The computed UTl time, on the other 
hand, is assumed to be consistent with the x axis of the implied TRF 
frame (see Section 4.3). 

The Celestial Reference Frame (CRF) employed for the integration of 
the equations of motion, is realized from the implied TRF frame through 
the following transformation (Mueller, 1969): 

(^) = (SNP)T(TRF) (2-30) 

with 

S = Ra (-xp) Rj (-yp) R 3 (GAST) (2-31) 

where according to the MERIT standards the following quantities have 
been used. 

P Precession matrix based on the lAU (1976) system of astronomical 
constants (Lieske, 1979) 

N Nutation matrix based on the 1980 lAU theory of nutation (Wahr, 
1981a). This matrix implies a pole whose nearly diurnal space-fixed 
and earth-fixed motions vanish. This pole is referred to as the 
Celestial Ephemeris Pole (CEP) (Mueller, 1981; Moritz and Mueller, 
1987) 

GAST ~ GMST (QhUTl) + f(UTl) + EQ.E (2-32) 

GMST(O^UTl) = 6 ^ 41"' 50?54841 + 8640184?812866 Tu + 0?093104 Tu^ 

- 6?2 X 10“« Tu® (2-33) 

where 

f = conversion factor from Universal time to sidereal time 

= 1.002737909350795 + 5.9006x10"“ Tu - 5.9x10"“ Tu® 

Tu = Julian centuries elapsed from J2000.0 

UTl = UTC(USNO) + (UTl - UTIR) + [UTIR - UTC(BIH)] + 

+ [UTC(BIH) - UTC(USNO)] (2-34) 
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UTC(USNO) = USNO Universal C<x>rdinated time (time scale used to 
time tag the observations of the GLTN stations) 

UTl - UTIR = Tidal variations of UTl caused by zonal tides with 
periods up to 35 days (BIH Annual Report 1981 onwards, Table Bl) 

UTIR - UTC(BIH) = Variations of the regularized UTl (i.e., UTIR) 
from the UTC(BIH) (BIH Annual Reports, Table 8) 

UTC(BIH) - UTC(USNO) =Variations of UTC(USNO) in relation to 
UTC(BIH) 

EQ.E = • cos (e + Ae) (2-35) 

= Nutation in longitude computed from the 1980 lAU nutation 
theory 

e = Mean obliquity of the ecliptic 

= 23* 26’ 21V448 - 46V8150Tu - 0'.’00059Tu* + 0V001813Tu3 (2-36) 

Ae = Nutation in obliquity computed from the 1980 lAU nutation 
theory 

The CEP pole positions Xp and yp in equation (2-31) have been taken 
from the smoothed values of Circular D (BIH Annual Report, 1983, 1984, 
Table 7). These pole positions are referenced to the 1979 BIH system 
during the first period of the MERIT Main Campaign (Sept. 1983 - Dec. 
1983), while during the remaining period of the campaign (Jan. 1984 - 
Dec. 1984) they are referenced to the BIH Terrestrial System (BTS). Our 
study in not affected by this transition because in shifting from the 
1979 BIH system to the BTS system a nonrotation constraint was applied 
to assure the continuity of the BIH system (Boucher and Feissel, 1984). 

2.2.5 Orbital Model 

The set of elements necessary to determine a satellite orbit constitutes 
the orbital model of the satellite. Thus, an orbital model consists of all 
those elements that are essential to formulate and integrate the 
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equations of motion (i.e., perturbations to be considered and associated 
assumptions, initial conditions, etc.). 

The Lageos state vectors required in the evaluation of the 
observation equations (Section 2.2.6) are obtained by numerically 
integrating the equations of motion (2-25). In this equation the inertial 
accelerations are expressed relative to the geocenter. The components, 
however, of Lageos* inertial acceleration caused by the nonsphericity of 
the earth are evaluated in the TRF frame and subsequently are 
transformed into the corresponding CRF frame, while the components of 
the remaining inertial accelerations are directly expressed in the 
corresponding CRF frame. The two-step procedure used to evaluate the 
Lageos inertial acceleration caused by the nonsphericity of the earth is 
necessary because the gravity potential coefficients are conveniently 
expressed in a TRF frame. Expressing these coefficients directly in a 
CRF frame would make them time dependent and therefore a potential 
source of unnecessary complications. 

The aim of the present study is not to estimate the Lageos orbit 
with the highest accuracy but rather to model it as simply as possible 
and yet be able to recover the baselines with an accuracy compatible to 
or even better than that of the observations. With this in mind the 
MERIT standards have been violated whenever the proposed model is 
complicated and cumbersome to incorporate into the solution. In such 
cases a simpler model has been adopted. It turns out, however, that in 
some cases the employed orbital model could be further simplified 
without affecting the accuracy of the recovered baselines (Section 4.6). 


2.2. 5.1 Point mass gravitational acceleration . The point mass Lageos 
gravitational acceleration based on the effects of the three major 
perturbing bodies (Earth (E), Moon (M) and Sun (S)) and expressed 
relative to the geocenter takes the following form (Cappellari et al., 
1976; Pavlis, 1982) 









(2-37) 


31 


where 


=: Lageos position vectors relative to the center of 

mass of the Earth, the Moon and the Sun respectively 

=: Earth position vectors relative to the center of mass 
of the Moon and the Sun respectively 

=: ratios of lunar and solar masses to the mass of the 
Earth 

=: geocentric gravitational constant 

For the evaluation of equation (2-37), one needs the Lageos geocentric 
I>osition vector as well as the geocentric position vectors of the Moon 
and the Sun respectively. The Lageos geocentric position vectors are 
obtained from the numerical integration of the equations of motion while 
the heliocentric position vector of the Earth and the geocentric position 
vector of the Moon are calculated from the information supplied by the 
DE/LB200 lunar planetary ephemeris (Standish, 1981). This ephemeris is 
disseminated in terms of Chebychev coefficients. These coefficients can 
only be used to calculate the geocentric positions of the Moon and the 
barycentric positions of remaining planets and the Sun. With this 
information, however, one can very easily calculate the position of any 
desired planet with respect to any of the remaining planets and to the 
sun as well. 

The reference frame implied from the computed coordinates of the 
planets has been accurately adjusted to the dynamical equinox J2000.0 
(ibid.) The Chebychev coefficients of the DE/LE200 lunar planetary 
ephemeris are based on the planetary coordinates estimated through the 
numerical integration process involved in the adjustment of 
interplanetary observations collected over a long period of time (ibid.). 
In this adjustment, the planetary equations of motion were formulated on 
the basis of the isotropic, parametrized post-Newtonian (PPN) n-body 
metric (Moyer, 1971). The independent variable in the PPN metric is the 
Barycentric Dynamical Time (i.e., coordinate time), and therefore this 
time scale should be used as an entry to the DE/LE200 ephemeris. 




fg. r| 
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The TDB time at any epoch j is computed from the ephemeris time of the 
same epoch via the following formula 

TDBj = tj + AT (2-38) 

where 

TDBj = Barycentric Dynamical Time at the epoch j 

tj = ephemeris time at epoch j, obtained from equation (2-24) 

In the PPN framework the ephemeris time coincides numerically but not 
conceptually with the Terrestrial Dynamical Time (i.e., proper time). 
ThuSf TDB time at any epoch is obtained by adding to the ephemeris 
time of the same epoch a small correction AT. This correction accounts 
for the general relativistic effects involved in the transformation of 
proper time (i.e., TDT time) to coordinate time (i.e., TDB time). An 
approximate value for the correction AT is given by the following 
formula (Astronomical Almanac, 1984) 

AT = 0?001658 sin (g) + 0.000014 sin (2g) (2-39) 

where 

g = 357:53 + 35999.05 Tu (2-40) 

In both of these equations (2-39 and 2-40), higher-order terms have 
been neglected, g designates the mean anomaly of the Earth in its orbit, 
and Tu designates the Julian centuries elapsed since J2000.0. 

To complete the evaluation of equation (2-37), one still needs the 
ratios of the lunar and solar masses to the mass of the earth as well as 
the geocentric gravitational constant. For the mass ratios, we have 
used the values recommended by the MERIT standards, but for the 
geocentric gravitational constant the value estimated simultaneously with 
the potential coefficients of the PGS1680 gravity field has been used: 

GMe = 3.986004359 x 10** mVs» (2-41) 
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The scale in the range dynamic mode methods is implied not only by the 
adopted value of the geocentric gravitational constant through the 
modified Kepler’s third law but also through the speed of light used to 
convert time measurements to range measurements. Thus, the adopted 
value of the geocentric gravitational constant should also be consistent 
with the speed of light implicit in the range observations. In the 
present study we have used the speed of light proposed by the MERIT 
standards (i.e., c = 299,792,458 m/s) (Lerch et al., 1985; Christodoulidis 
et al., 1985). 

The associated partial derivatives of equation (2-37) contributing to 
the variational equations of state (i.e., to matrix Y(t), equation (2-28)) 
are given in (Cappellari et al., 1976, eq. 4-21; Pavlis 1982, eq. 13). 

2.2.5.2 Gravitational acceleration due to nonsphericity of the 
gravitational potential . The inertial acceleration induced on the satellite 
by the nonsphericity of the earth is obtained via the gradient of the 
perturbing potential. The perturbing potential is a scalar function 
describing the nonspherical part of the geopotential in terms of an 

infinite spherical harmonic series (Heiskanen and Moritz, 1967): 

GMcf • fapl” n 

VMs(r.^.^) = - 7 - I l~J I [Cnm cos(inX) + S„m sin(mX)]Pn„(sin4>) 

” ln=2 ^ m=o J 

(2-42) 

The zero-degree harmonic has been modeled in equation (2-37) and 
therefore is not included in the above equation. The first-degree 
harmonics are also not included because the origin of the PGS1680 

implied (TRP) coincides with the the center of mass of the earth (see 
Section 2.2.4). With the gradient of the perturbing potential the 
components of Lageos’ inertial acceleration, caused by the nonsphericity 
of the earth, are expressed in the PGS1680 implied (TRF) frame. 

Incorporation of this acceleration into the equations of motion (2-25) 

requires transformation of its comiK>nents from the (TRP) frame into the 
corresponding (CRP) frame via the transformation equation (2-30) 
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Rns = (SNP)T r 


(2-43) 


where 
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(2-44) 


is the gradient of the perturbing potential function V^s (equation 2-42). 
The expressions for the partial derivatives of the perturbing potential 
function are given in (Cappellari et al., 1976). The radius (ag) of 
the reference sphere, also needed in the evaluation of the perturbing 
I>otential, is the same with the radius employed in the estimation of the 
PGS1680 gravity field (i.e., a^ = 6378144.11 m). For our study we have 
truncated this field at degree and order 12 because perturbations 
caused by higher harmonics over a two-week period contaminate the 
computed SRD observables with errors having magnitude well below the 
noise level of the SRD quasi-observables (see Section 4.6). Furthermore, 
a nonvariant nature of the coefficients C 21 and Sa, has been adopted, 
although it is well known that these two coefficients are largely affected 
by the forced diurnal motion of the figure axis caused by the 
nonrigidity of the earth (Moritz and Mueller, 1987). 

The associated partial derivatives of equation (2-42) contributing to 
the variational equations of state are given in (Cappellari et al., 1976, 
eq. 4-54; Pavlis, 1982, eq. 21). 


2.2.S.3 Lageos tidal inertial acceleration . According to the MERIT 
standards, the tidally induced variations in the earth's external potential 
should be incorporated in the orbital model as variations in the 
geopotential coefficients (Melbourne et al., 1983; Banes et al., 1983). In 
order to save computing time a two-step procedure is proposed to carry 
out the implementation of these variations (ibid). In the first step the 
variations of the geopotential coefficients are computed on the basis of a 
nominal frequency independent Love number k„, while in the second 
step these variations are corrected to account for the frequency 
dependent nature of the nominal Love number k„. 
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The frequency dependent variations of the Love number k„ have been 
estimated for an elliptical, rotating, elastic fluid outer core and solid 
inner core, oceanless earth (Wahr, 1981b). Although the proposed 
two-step procedure is computationally more effective than any one-step 
procedure, it is not appropriate for our investigation because one still 
needs to evaluate many trigonometric functions at each of the observing 
epochs (Melbourne et al., 1983). This however, not only would make the 
orbital model more complicated but also it would make the 3RD method 
computationally less efficient. Therefore, it was decided to compute the 
tidally induced space potential by assuming a solid earth (i.e., oceans 
not included) which exhibits the same elastic response over all possible 
orders within a certain degree (Diamante et al., 1972; Pavlis, 1982). With 
such an earth model the tidally induced potential on the surface of the 
earth takes the following form (Diamante et al., 1972; Goad, 1977; Pavlis, 
1982): 

Ut„ = Z kn Ut„ (se) (2-45) 

D „=2 n 

where k„ is the nominal Love number of degree n and Uy (se) is the n^^ 

n 

surface harmonic of the tidal potential. Solving the Dirichlet problem 
the tidally induced space potential is obtedned: 


Ut 


D 


CD 

z 

n = 2 



kn Ut^Cue) 


(2-46) 


where |r| is the norm of the Lageos position vector expressed relative 
to the center of mass of the earth. To the first order, terms with n > 2 
in equation (2-46) can be neglected (Diamante et al., 1972), and therefore 
this equation takes the following form 


Ut 


D 


a. 

-A ka Ut^ (ue) 

R 


(2-47) 


where 
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(2-48) 


The quantity k 2 denotes the second-degree Love number while Mb 
denotes either the lunar or the solar mass> or for that matter the mass 
of any planet that is considered to be a disturbing body. In this study 
the Moon and the Sun have been considered as the only disturbing 
bodies. The vectors R and Rb designate the geocentric position vectors 
of the satellite and the disturbing body respectively. The components 
of these vectors are expressed in the corresponding (CRP) frame. In 
this frame the tidally induced acceleration on Lageos takes the following 
form (Diamante et al., 1972; Pavlis, 1982): 


- 3 GMb r / ^ 1 

Rtd = 2 p' j; ■ TT77 [(l - 5(Ub • u)*)u + 2(ub • u)Ubj (2-49) 

|Rbl |R| 

where 

Ub = and u = -^ (2-50) 

iRbl IrI 

To account for a phase lag produced by the earth’s dissipative forces, 
the vector Rb in equation (2-50) has been replaced by another vector 
R*. This vector is obtained from the vector Rb via the following 
transformation 

RJ - R3(“fiL)^b (2-51) 

where (=0*.35) is the phase lag. The R 3 rotation is performed about 
the third axis of the corresponding CRP frame (see Section 2.2.4). In 
equations (2-47) and (2-49) the value 0.29 was adopted for the 
second-degree Love number ka. This value is different from that 
proposed by the MERIT standards (ka = 0.30). This deviation, although 
not of much importance, is justified since the tidal corrections applied 
in the estimation of the PGS1680 gravity field were based on the altered 
value of k 2 (i.e., k 2 = 0.29). The permanent tidal deformation affecting 
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the potential coefficient C 20 is inherently present in equation (2-49). 
Consistent incorporation of this equation in the equations of motion 
requires that the permanent tidal deformation is not included in the 
PGS1680 C 20 value. This, however, seems to be the case for the 

PGS1680 gravity field (Melbourne et al., 1983; Christodoulidis et al., 
1985). Furthermore, the ocean tidal perturbations are not included in the 
orbital model of Lageos not only because they are small (i.e., one order 
of magnitude smaller than the solid earth tidal effects (Section 4.6)) but 


also because their evaluation 
computations considerably. 

would 

increase 

the bulk 

of 

the 

The contribution of the 

tidally 

induced 

acceleration 

to 

the 

variational equations of state is 

given in 

(Pavlis, 1982, eq. 30). 

In 

that 

equation the vector Rb should 
equation (2-51). 

be replaced by 

the vector 

Rg 

from 


2.2.S.4 Lageos solar radiation pressure acceleration . The 
acceleration induced on Lageos due to photon momentum transfer is 
referred to as solar radiation pressure acceleration, and it is given by 
the following formula (El'Yasberg, 1967; Cappellari et al., 1976; Pavlis, 
1982). 

Rsr - r • c 

The eclipse factor 7 assumes the values zero, one, or any other value in 
between depending on whether the satellite is in complete shadow 
(umbra), in sunlight, or in partial shadow (penumbra) respectively. In 
our study the eclipse factor 7 is determined by a simple cylindrical 
model (Cappellari et al., 1976; Pavlis, 1982). This model is easy to 
incorjxjrate into the equations of motion, but it does not differentiate 
between umbra and penumbra regions. A full model for the earth's 
shadow, as proposed by the MERIT standards, would increase the bulk 
of computations thereby complicating the solution. It is rather doubtful 
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(2-52) 
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that this complication would make any difference. The mean solar flux S 
is the amount of photon energy flow through a unit surface per unit 
time at a distance of one astronomical unit (AU) (i.e., AU r 1.4959787066 
X 10** meters) (Melbourne et al., 1983). The ratio (S/C) is the photon 
momentum transfer to a unit surface per unit time at a distance of one 
astronomical unit. The value (4.5605 x 10“* N/M*) was used in our 
study for the ratio (S/C) as proposed by the MERIT standards. The 
position vectors R and Rs in equation (2-52) designate the geocentric 
position vectors of Lageos and the Sun respectively. The reflectivity 
coefficient (Cr) depends not only on the mechanism of light reflection 
but also on the thermal emission distribution of the satellite surface. 
The monthly values for the reflectivity coefficient (Cr) are shown in 
Table 1 for the entire MERIT campaign. These values have been 
estimated together with other parameters in the adjustment of the MERIT 
laser range data performed by the GEODYN II programs (Pavlis, 1986, 
private communication). In the present study we have used the 
reflectivity coefficient values listed in Table 1 instead of using the 
value proposed by the MERIT standards. 


Table 1 Lageos Along— Track Acceleration and Its Reflectivity 
Coefficients 



Magnitude of Lageos Along-Track 
Acceleration x 10~*^ m/s* 

Lageos Reflectivity 
Coefficients 

Sep. 

1983 

-2.909 

1.141 

Oct. 

ft 

-3.549 

1.136 

Nov. 

tt 

-3.893 

1.135 

Dec. 

ft 

-3.825 

1.133 

Jan. 

1984 

-4.343 

1.136 

Feb. 

tt 

-4.319 

1.134 

Mar. 

ft 

-4.065 

1.109 

Apr. 

tt 

-3.550 

1.057 

May 

tt 

-3.189 

1.096 

Jun. 

tt 

-3.393 

1.139 

Jul. 

tt 

-3.928 

1.132 

Aug. 

tt 

-2.946 

1.126 

Sep. 

tt 

-2.301 

1.126 

Oct. 

tt 

-2.524 

1.126 
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This deviation of the MERIT standards, although plausible, will not affect 
the accuracy of the estimated baselines. The effective area (A) of the 
surface normal to the incident light is given, for a spherical satellite 
like Lageos, from (El’Yasberg, 1967) 

A = ttRs’ (2-53) 

where Rg is Lageos’ radius (i.e., Rj = 0.30 m) (Melbourne et al., 1983). 
Finally, the Lageos mass (M) of 407 kg has been used in the evaluation 
of equation (2-52). 

The associated partial derivatives of the solar radiation pressure 
acceleration contributing to the variational equations of state are given 
in (Cappellari et al., 1976, eq. 4-161 and 4-162; Pavlis, 1982, eq. 25 and 
26). 

According to the MERIT standards the inertial acceleration induced 
on the satellite due to the diffused reradiated light from the earth 
(earth albedo) is not included in the orbital model of the Lageos 
satellite. 

2.2.5.5 Lageos along-track empirical acceleration . Ever since the 
launch of the Lageos satellite it has been observed that its semimajor 
Gods decreases at a rate of 1 mm/day. This has been traced to an 
unexpected and still physically unmodelled along-track acceleration 
acting on the Lageos satellite. Attempts to explain the origin of this 
mysterious acceleration have either totally or partially failed. These 
attempts are based on a variety of possible causes ranging from 
assuming helium concentrations at satellite altitudes (Rubincam, 1980) to 
considering the solar eclipses (Rubincam et al., 1985). Although all of 
these attempts have partially failed, it is quite clear that this 
acceleration is the result of a combined effect caused by the 
asymmetries of the earth’s albedo and by the charged particles traveling 
in the vicinity of Lageos (Smith et al., 1985; Alfonso et al., 1985). Since 
the physical process producing this acceleration is unknown, its 
modeling has been accomplished with an empirical model. This model 
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assumes a "drag-like" force acting against the satellite motion (Pavlis, 
1982): 



(2-54) 


where « is the magnitude of the along-track acceleration. The monthly 
magnitudes of this acceleration are also listed in Table 1 for the entire 
MERIT campaign. These values have been estimated with the GEODYN II 
program (Pavlis, 1982, private communication). 

Contributions to the variational equations of state due to this 
acceleration are neglected because of their small magnitudes. 


2.2.6 Normal Equations 

The observation equations for the SRD quasi-observables are readily 
obtained through a Taylor series expansion of equation (2-22). In this 
expansion only the zero- and first-order terms are retained while all of 
the remaining higher-order terms are neglected. The expansion is 
performed about the approximate earth-fixed station coordinates and the 
celestial initial state vector of the corresponding arc: 




(2-55) 
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- (dYi , dY^ , dRo , dRo);xi: 


(2-59) 


Lj = the computed minus the generated SRD quasi-observable 

Vj = residual corresponding to the j^h sRp observable 

The matrices (Bj)ix6 and (Sj)ix 3 are readily obtained by differentiating 
equation (2-22), while the state transition matrix [Yj(t)] 3 x 6 is obtained 
in the numerical integration process of the variational equations of 
state. The celestial satellite coordinates at the epoch j needed to 
evaluate the vectors (Bj)ix6 and (Sj)ix 3 as well as the scalar Lj are 
obtained from the numerical integration of Lageos* equations of motion 
(see Sections 2.2,3 through 2.2.5). The adjusted parameter vector Xj 
contains corrections to the earth-fixed approximate coordinates of 
stations 1 and 2 (i*e., dYi and dY 2 ) and to the corresponding celestial 
initial state vector (i.e., dR© and dRo). Extension of equations 
(2-55)-(2-59) to include all the available SRD observables and all the 
observed satellite arcs leads to the following equations 


V = A*X + L 
where 


(2-60) 


A* = [ B* ; C* ] 


(2-61) 


B* 


C* = 



'B6p 




JY . 

nx 3 k 



B6p 


3R 


3R 

nx 3i 



(2-62) 


(S)„x3i ■ [Y*(t)]3^,,„ (2-63) 


^ 3ix 


X = [dY , d(lR^,Isj] 


(2-64) 


42 



L - vector containing the computed minus the SRD quasi-observables 
V = residual vector 

Vector fip contains all the available SRD observables, vector Y the 
Cartesian coordinates of all the observing stations, vector Z the 
Cartesian coordinates of the observed satellite positions, and vectors 
Ilfo» the initial state vectors of all the arcs involved. The adjusted 

A 

parameter vector x contains the corrections to the approximate 
earth-fixed coordinates of the observing stations (i.e., dY) together with 
the corrections to the initial state vectors of the observed satellite arcs 
(i.e., The integers i and n denote the number of the 

observed satellite positions and the number of observations, while the 
integers k and m denote the number of the observing stations and the 
number of the observed satellite arcs respectively. 

A close examination of equations (2-62) and (2-63) reveals that the 
submatrices (B„x 3 k) and (Snxsi) would be exactly the same even if the 
satellite positions were treated as auxiliary independent points in space 
(i.e., geometric approach). The constraints imposed on the observed 
satellite positions to lie in the corresponding satellite arcs are applied 
through the state transition matrix (Y*(t) 3 ixsm)* 

Singularity A (see Section 2.1.3) affects the dynamic and geometric 
solutions in exactly the same way because the submatrix B* in equation 
(2-61) is the same for both the geometric and the dynamic approach. 
With SRD observables singularity A occurs not only from the resulting 
geometry of one station and its observed targets (see Section 2.1.3) but 
also from the geometry of two coobserving stations and their observed 
targets. Singularity B or singularity C cannot exist in a dynamic 
solution because the structure of the matrix (S„x 3 i) is altered by its 
multiplication with the state transition matrix (Y*(t) 3 ix 6 m)* The 
alteration of the matrix (S„x 3 i) not only differentiates the dynamic from 
the geometric approach but also furnishes the dynamic approach with 
better stability characteristics (see Sections 4.2 and 4.4). 

Taking into consideration that the state transition matrix is different 
from epoch to epoch, one can readily prove that in the absence of 
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singularity A the design matrix A* (apart from the ill-defined origin of 
longitudes) is nonsingular if the adjusted parameter vector consists only 
of corrections to the approximate station coordinates and to the initial 
state vectors. This is not surprising because the TRF frame, with an 
ill-defined origin of longitudes, is implied by the PGS1680 gravity field 
while the CRF frame is subsequently realized via the transformation 
equation (2-30). Including polar motion and/or variations of UTl in the 
adjusted parameter vector results in an extremely ill-conditioned design 
matrix A* because polar motion and station coordinate are nearly 
inseparable, while variations in UTl and in the satellite node are 
inseparable parameters as well (Van Gelder, 1978; Pavlis, 1982). 

In the present study, polar motion and the variations in UTl are 
not included in the adjusted parameter vector. Thus, after resolving 
the problem of the ill-defined origin of longitudes (see Section 4.4) we 
proceed with the formation and the solution of the normal eqiiations. 
Using the same arguments as in Section 2.1.2 and the observation 
equations (2-60), the normal equations take the following form (Uotila, 
1987) 

(A*TPA* + Px)X + A*TPL = 0 (2-65) 

where P and Px are the weight matrices associated with the SRD 
observables and the adjusted parameter vector respectively. Since the 
weight matrix P is diagonal the normals are formed sequentially through 
the following formula 


n ATA. 

A*TpA* + Px = I + Px 


( 2 - 66 ) 


where Aj (i.e., equation (2-56)) is the j^^ row of the design matrix A*, n 
is the total number of the SRD observables and <rj is the variance for 
the j^^^ SRD observable. These variances are computed via the following 
formula 
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where and are the variances of the actually observed and the 

interpolated ranges respectively. The inversion of the normal equation 
matrix (2-66) was obtained with the Cholesky algorithm (Uotila, 1967). 



Chapter 3 

GENERATION OF THE OBSERVABLES 


This chapter starts with a description of the SLR system in an attempt 
to identify and understand the systematic errors affecting the laser 
ranges. It continues with a description of the data set employed in this 
investigation and finally ends with the generation of the simultaneous 
range and SRD observables. These two observables constitute the input 
to the geometric and SRD methods respectively. 


3,1 SATELLITE LASER RANGING 

A satellite laser ranging system consists of three basic components: 

(i) the ground segment, 

(ii) the atmospheric channel, and 

(iii) the spaceborne segment. 

The ground segment consists of a global network of fixed and 
highly mobile satellite laser ranging stations forming a network 

configured to allow measurements of the plate tectonic motions (Coates et 
al., 1985). Tectonic plate motions are essential in understanding the 
geodynamic processes necessary for earthquake and volcano erruption 
predictions. Each of the stations in the network is equipped with the 
necessary hardware to produce, emit, receive and measure the round- 
trip flight time of very short laser pulses to a retroreflector equipped 
artificial satellite such as LAGEOS. 

The atmospheric channel is the optical path followed by a laser 
pulse in its round trip from the station to the satellite. 

The spaceborne segment consists of approximately 14 retroreflector 
equipped satellites (Degnan, 1985). For geodesy and geodynamics Lageos 
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is an example of such a satellite in orbit high enough not to be 
influenced by the difficult to model high frequency variations of the 
gravity field and the atmospheric drag but yet at low enough altitude to 
assure good signal returns to the tracking stations. Therefore, the 
propagation of the orbital errors in the estimated geodetic parameters is 
substantially reduced. This error reduction is very important because 
variations in certain geodetic parameters such as baselines, polar motion, 
and length of day are routinely used in understanding the mechanisms 
driving geodynamic processes. 

In the operational environment, depending on the technology 
employed and the models used, each component of the satellite laser 
ranging system will contribute in part to the total error affecting the 
inferred geometric range. The next section contains, for each component 
of the SLR system, a brief discussion of its operational principles, the 
error sources, their status during the MERIT Main Campaign and the 
future possibility of either reducing or eliminating them. 


3.2 SATELLITE LASER RANGING SYSTEM, ITS COMPONENTS AND THEIR 
CONTRIBUTION TO THE TOTAL ERROR BUDGET 

3.2.1 Hardware of the Ground Segment 

For each satellite ranging system, the hardware of the ground segment 
consists of the laser transmitter, the laser receiver, their transmittimg 
and receiving optics, the timing subsystem and the computer. 

The laser transmitter in most of the modern laser systems consists 
of a mode-locked Nd:YAG laser oscillator followed by one or more Nd;YAG 
laser amplifiers. The name Nd:YAG is derived from the crystal used in 
the light amplification by stimulated emission of radiation (i.e., lasing 
process) which is a YAG crystal (Yttrium aluminum garnet : Y 3 AI 5 O 2 ) 
"doped" with Neodymium (Johnson et al., 1978). Since the mode-locked 
Nd:YAG lasers operate in a single spatial mode they are not affected by 
"wavefront-distortion" errors (Degnan, 1985). The biases introduced by 
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the mode-locked transmitters are of the order of subcentimeter level. 
The crystals used in the lasing process, for the stations participating in 
the MERIT Main Campaign, are reported in the SLR coordinator’s report 
and its updates (Schutz, 1983b). 

The laser receiver is designed to measure the round-trip flight time 
of the laser pulse to a retroreflector equipped satellite. This time 
interval is multiplied by the speed of light and divided by two to infer 
the optical range from the station to the satellite. The basic elements of 
a laser receiver are the photomultiplier, the discriminator and the time 
interval unit. 

The photomultiplier is a device used to detect the incoming laser 
pulse. Its principle of operation is based on the photoelectric effect 
(Halliday et al., 1962; Drain, 1980). Most of the SLR systems 
participating in the MERIT Main Campaign made use of the conventional 
type photomultipliers referred to as dynode-chain photomultipliers 
(Degnan, 1985). The time it takes for the photoelectrons to propagate 
from the photocathode to the anode via the dynodes is called transit 
time. If the transit time were constant it could be completely accounted 
for through either calibration or common channel procedures (ibid). 
Variations in the transit time, referred to as transit time jitter, 
influence the inferred ranges by as much as 15 cm (ibid). This error is 
mainly caused by the motion of the satellite image within the 

photocathode when the instrument tracks the satellite. However, 
successful focusing of the satellite image onto the photocathode reduces 

this error to the 1 cm level. Other factors such as the impulse 

response of the PMT’s, the amplitude of the input signals and the 
background radiation also contribute to this error. These problems are 
currently being solved with the replacement of the conventional PMT’s 
with the so-called microchannel plate photomultiplier tubes (MCP/PMT) 
recently introduced on the market. These photomultipliers are 

characterized by well-defined photoelectron path lengths with much 
shorter transit times, and greatly reduced sensitivity in the image 
position effects, the strength of the input signal and the background 
radiation. When using the MCP/PMT photomultipliers, the resulting 
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errors in the inferred range appear to be below the 5 mm level (ibid). 

The purpose of a discriminator is to define on the photomultiplier’s 
output waveform a timing point and subsequently to generate a 
rectangular logic pulse that starts or stops the time interval unit. The 
output waveform has a quasi-Gaussian form with a randomly varying 
amplitude. This amplitude variation introduces, in the determination of 
the timing point, a time bias which is highly repeatable and can be 
estimated if the amplitude of the input pulse is measured and recorded 
along with each observation. In practice the amplitude-dependent time 
bias is determined experimentally and is compensated for by 
incorporating a hardwired circuitry into the discriminator. The degree 
of success in implementing this circuitry is determined experimentally 
and is shown in the time walk characteristic of the discriminator. The 
time walk characteristic is a curve obtained by plotting the signal 
amplitude dependent time biases versus the input signal amplitudes. 
The time walk characteristic shows the time bias introduced by 
amplitude variations, while its RMS deviation from the zero horizontal 
line characterizes the efficiency of the discriminator (ibid). For the 
discrimators used during the MERIT Main Campaign, this RMS deviation 
reached the value of 1.5 cm, while for discrimators currently appearing 
on the market this value has been reduced to the 0.2 cm level. The 
latter discriminators are currently being tested for implementation in the 
continuously upgraded SLR systems. 

The purpose of the Time Interval Unit (TIU) is to measure the 
round-trip flight time of the laser pulse. The rectangular logic pulse 
generated by the start discriminator activates the time interval counter 
while the corresponding logic pulse from the stop discriminator 
commands the counter to stop. The basic component of the TIU is an 
oscillator which determines the stability and accuracy of the TIU. The 
oscillators used by the SLR stations which participated in the MERIT 
Main Campaign were either cesium beam type or rubidium type (Schutz, 
1983b). To achieve maximum accuracy, the measurement of the round-trip 
flight time is split up in three parts (i.e., T = T1 + T12 - T2), where T1 
is the time elapsed from the starting epoch to the first following 
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positive crossover of the oscillator, T2 is the time elapsed from the 
ending epoch to the next following positive crossover, and T12 is the 
time interval between the eiforementioned positive crossovers. T12 is 
obtained by multiplying the number of intervening positive crossovers 
(N12) by the period (TO) of the master oscillator (i.e., T12 = N12 x TO). 
The fractional times T1 and T2 are accurately measured either by 
charging and discharging a capacitor with constant but different 
currents or by using a second oscillator which is slightly off from the 
master oscillator. Common biases introduced in measuring T1 and T2 are 
canceled out because the times T1 and T2 are subtracted in the 
computation of the round trip flight time (i.e., T). Residual errors at 
the cm level are still present. These errors can be reduced at the mm 
level with the use of streak-cameras employed in the Optical Time 
Interval Unit (OTIU) currently being investigated for implementation in 
the new two-color laser receivers (Abshire et al, 1985). 

The transmitting optica are used to align the laser pulse towards 
the satellite being tracked, while the receiving optics are used to 
receive and focus the reflected laser pulse onto the cathode of the 
photomultiplier. Unsuccessful focusing introduces the image position 
effects previously mentioned. A substantially reduced single 
dual-purpose telescope performs both of the above functions, thanks to 
the technological advances in the field of signal detectors, 
photomultipliers and discriminators. These advances introduced a 
substantial reduction in the station design which in turn triggered the 
construction of the highly transportable laser ranging systems TLRS-1 
(Silverberg, 1982; Shelus, 1983), TLRS-2 (Transportable Laser Ranging 
System No. 2), etc. These systems are extremely valuable in the study 
of geophysical processes because of their ability to make observations of 
limited time in remote areas, in a hostile environment. 

The time receiver is either a LORAN— C or a GPS receiver and is 
used to time-tag each observation in the Universal Coordinated Time 
(UTC) scale. The synchronization accuracy with a LORAN-C receiver is 
of the order of 1 ^ts, while with a GPS receiver this accuracy is of the 
order of 50 ns. Synchronization errors affect the inferred optical range 
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by about 4 mm//is. The laser stations used in this study belong either 
to the Goddard Laser Tracking Network (GLTN) or to the Participating 
Laser Network (PLN) (Shawe et al., 1985), and their observations are 
time-tagged with the UTC scale kept by USNO. According to the MERIT 
standards and for reasons explained in Chapter 2, the UTC(USNO) has 
been transformed to UTC(BIH). 

The computer is used as an auxiliary equipment to control satellite 
tracking operations, to assist the operator with such functions as data 
quality and quantity assessment, maintenance, testing procedures, etc. 

3.2.2 Atmospheric Channel 

As the laser pulse propagates through the atmospheric channel it 
experiences a continuously varying refractive index. This variation 
depends primarily on the variations of the local pressure with only a 
weak dependence on the local temperature and local humidity (Degnan, 
1985; Abshire, 1985). A varying refractive index, on the other hand, 
bends the laser pulse according to Snell’s law and also decreases the 
group velocity of the laser pulse as it travels through lower pressure 
layers at higher altitudes. The error due to bending of the laser pulse 
is relatively small and reaches a maximum value of 3-4 cm at 10 degrees 
elevation while the error due to the decrease of the group velocity is 
very large, reaching the value of about 13 m at the same elevation 
(Abshire, 1985). A great number of formulas have been developed to 
correct the inferred optical length of the laser pulses (i.e., the inferred 
laser ranges) for atmospheric refraction effects. In the present study 
and according to the MERIT standards, the Marini and Murray formula 
has been used to correct for these effects (Marini and Murray, 1973). 
This formula is based on the assumption of a spherically symmetric 
atmospheric refraction and it uses only the pressure, temperature and 
relative humidity taken at the ranging site. This formula is in error at 
the 4-6 cm level as the satellite reaches an elevation angle of 20 
degrees (ibid). Today, use of two-color laser ranging systems equipped 
with streak-camera receivers promises an atmospheric refraction 
correction with an accuracy down to the mm level, thanks to the weak 
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dependency of the group refractivity on the water vapor at optical 
wavelengths (Abshire, 1985). 

3.2.3 Space Segment 

Although the space segment of the SLR system consists of several 
satellites, the LAser GEOdynamics Satellite (LAGEOS) is devoted 
exclusively to geodynamic and geodetic applications (Moritz and Mueller, 
1987). Lageos is a passive sphere with a 59.988 cm diameter orbiting 
the earth at an altitude of about 5900 km. Its mass-to-area ratio of 
1.44x10® kg/m® effectively minimizes the solar radiation pressure and 
atmospheric drag perturbations. The high altitude of Lageos’ orbit not 
only reduces the effects of the poorly modeled high frequency 
variations of the gravity field but also warranties good simultaneous 
tracking of continental extent. The altitude of the orbit, however, is 
low enough to assure the geometric strength necessary for successful 
implementation of simultaneous laser satellite tracking methods. 
Consequently, only laser range observations to Lageos were facilitated in 
order to investigate the effectiveness of the 3RD and geometric methods 
in baseline determinations. 

The surface of the Lageos satellite is speckled with 422 

solid-cube-corner reflectors (OCR’s) made of fused silica and four made 
of germanium (Cohen et al., 1985). When the direction of the incoming 
laser beam relative to the normal of each individual CCR reaches the 
value of 25 degrees, reflection ceases to take place (Degnan, 1985). As 
a result, 10 to 15 OCR’s contribute to the laser pulse detected at the 
receiver. Therefore, it is difficult to locate for each returning pulse its 
reflection point which constitutes the ending i>oint of the inferred 
optical laser range. The location of the reflection point is needed to 
compute the correction necessary to transform the ending point of the 
optical laser range to the center of mass of the satellite. This 
correction is referred to as center of mass correction, and its value has 
been determined experimentally for different pulse widths prior to the 
Lageos launch. The standard deviation in estimating this correction is 
about 2 mm. In the current investigation the value of 24 cm has been 
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adopted for the center of mass correction. This value is supplied on 
each data record in the tape containing the observations (see Section 
3.4). Furthermore, the interference of the individual OCR returns, at 
the receiver’s level, may introduce a random error in the inferred laser 
range the standard deviation of which reaches the value of 1.15 cm 
(Pitzmaurice et al., 1977). This error is referred to as the coherent 
fading effect. 

3.2.4 Instrument Origin 

Effective use of laser range observations necessitates a clear 
identification of the starting and ending points of the inferred ranges. 
As already mentioned, the ending point is identified with the center of 
mass of the satellite. This is a natural choice because the equations of 
motion of the satellite are conveniently expressed relative to this point. 
The starting point is identified with a fixed reference point within the 
laser instrument and is referred to as instrument origin. The 
instrument origin usually coincides with the intersection of the 
telescope’s azimuth and elevation axes, but other points within the laser 
instrument may be used as well. Realization of the instrument origin is 
achieved either through calibration or through the common channel 
receiver approach (Abshire et al., 1984; Degnan, 1985). During the 
MERIT Main Campaign, calibration procedures were employed to identify 
the instrument origin through the estimation of the system delay 
(Schutz, 1983a). The system delay is measured by making repeated 
observations to a calibration target of known distance, usually before 
and after each satellite pass. With this information, the system delay 
introduced by the instrument’s electronics can be readily estimated. 
The distance to the calibration target is measured with a geodimeter 
located very close to the laser ranging instrument. Additional surveys, 
therefore, are necessary to determine the position of the geodimeter 
relative to the instrument’s origin. In this process, errors of the 
order of 2 cm may be introduced. In order to reduce these errors, 
some laser instruments are equipped with fiberoptics allowing for 
self-calibration (Silverberg, 1982). The common channel-receiver 
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approach, on the other hand, eliminates the need for calibration since 
the electronic system delay, except for a calibratable signal-amplitude 
effect, cancels itself out (Abshire et al., 1984; Degnan, 1985). This 
approach is currently being tested for implementation in the laser 
ranging systems. 

From the above discussion it is obvious that on the basis of the 
technology employed and the models used, one could come up with a 
standard deviation depicting the accuracy of the laser range 
observations recorded by a certain station. This approach, however, 
would not take into consideration errors resulting from improper 
calibration, from operator errors or from any other errors not being 
accounted for. Bearing this in mind, it was considered appropriate to 
estimate, for every station used in the present study, a standard 
deviation that would reflect the station’s overall performance during the 
MERIT Main Campaign. Such an estimate can be obtained by taking an 
average value of the monthly precision estimates determined for every 
station and for the entire MERIT Main Campaign by the University of 
Texas (Analysis of Lageos Laser Range Data, Sept. 1983-Oct. 1984). For 
the stations involved in our study, these estimates are shown in Table 2 
along with the station ID’s and the kind of laser instruments with which 
they were equipped during the MERIT Main Campaign. 
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Table 2. Station Location, Laser Instnnnents 
and Precision Estimates. 


NAME 

ID 

LOCATION 

LASER 

INSTRUMENT 

OBSERVATIONAL 
PRECISION (m) 

QUINC2 

7109 

Quincy, CA 

MOBLAS-8 

0.028 

MNPEAK 


Mount Laguna, CA 

MOBLAS-4 

0.033 

MAZTLN 

7122 

Mazatlan, Mexico 

MOBLAS-6 

0.12/0.05* 

GRF105 

7105 

Greenbelt, MD 

MOBLAS-7 

0.034 

PLATVL 

7112 

Platteville, CO 

MOBLAS-2 

0.125 

MCDON 

7086 

McDonald Obs., Ft. Davis, TX 

MLRS 

0.084 

TL0126 

7265 

Barstow, CA 

TLRS-1 

0.080 

OTAY 


Otay Mt. , San Diego, CA 

TLRS-2 

0.060 

QUINC3 

7886 

Quincy, CA 

TLRS-1 

0.070 

M0NPK2 


Mt. Laguna, CA 

TLRS-1 

0.060 

HOLLAS 

K 4 3 

Lure Obs . , Maui , HI 

HOLLAS 

0.042 

HUAHIN 


Huahine, Society Is., Pol. 

MOBLAS-1 

0.094 

ARELAS 

7907 

Arequipa, Peru 

AREfixed 

0.145 


^before and after upgrading. 


3.3 SYSTEMATIC CORRECTIONS OP THE OBSERVATIONS EXTERNAL TO 
THE SLR SYSTEM 

Effective application of the least-squares adjustment assumes constant 
adjusted parameters over at least the time span of the observations. 
Thus, baseline estimation from station coordinates requires corrected 
station coordinates for their temporal variations. Besides shocks and 
regional deformations, the temporal variations of station coordinates are 
caused either by tectonic plate motions or by earth tides. The regional 
deformations are ignored because either they are unknown or their 
effects are very small. For instance, the ocean loading effects for the 
stations used in this investigation are very small (Melbourne et al., 
1983). Since the time span of the observations covers about one year, 
the plate tectonic motions have also been ignored. Thus we have 
considered only the temporal variations of station coordinates caused by 
the earth tides. The tidal corrections are accounted for by correcting 
either the observations or the station coordinates. In the present study 
the traditional way of correcting the station coordinates has been 
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adopted. This correction has been conveniently formulated through the 
station displacement vector caused by the tidal deformation (ibid.). This 
formulation is based on the same elastic response of a solid earth over 
all orders within a certain degree, and if only the second degree is 
considered it takes the following form: 


Af* 



GM^r* 

GMeRJ^ 


[3A*(Rf-?)]Rj + 




where 



GMg 

f,r 

h2 

^2 


RsC-^l) • Rj (3-2) 

phase lag caused by the earth* s dissipative forces 
gravitational parameter of the attracting body. In the 
present study only the moon (j =2) and the sun (j =3) 
have been considered 
geocentric gravitational constant 

unit vector and the magnitude of the geocentric vector 
of the moon (j =2) and the sun (j =3) respectively 
unit vector and the magnitude of the stations* geocentric 
vector 

nominal second-degree Love number 
nominal second-degree Shi da number 

unit vector and the magnitude for the geocentric vector of 
of the moon (j =2) and the sun (j = 3) in the absence 
of dissipative forces 


The way the tidal displacements have been incorporated in our study 
differs in two aspects from what was suggested by the MERIT 
standards. The phase lag caused by the dissipative forces has been 
modelled in equation (3-2), while the second-degree Love and Shida 
numbers have been assumed to be frequency independent. The latter 
assumption results in a maximum error of 1.3 cm in the stations* height 
(ibid.). This assumption is well justified not only because the resulting 
error is well below the noise level of the observations but also because 
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it reduces the bulk of the computations considerably. 


3.4 DESCRIPTION OP THE DATA SET UTILIZED IN THIS INVESTIGATION 
Both the SRD and the geometric methods require strict simultaneity. 
Although baseline estimations based on exclusively simultaneous 
observations are largely insensitive to the orbital errors (Christodoulidis 
et al., 1981) and to reference frame model errors (Pavlis and Mueller, 
1983), no specific campaign was ever devoted to coordinate simultaneous 
tracking. This, together with the inability of the SLR systems to track 
a satellite through a cloudy atmosphere, makes it even more difficult to 
achieve extensive simultaneous laser tracking. 

Fortunately enough, as early as 1978 the lAU Symposium No. 82 on 
Time and the Earth’s Rotation" recommended setting up a working 
group to organize a program of international collaboration to Monitor the 
Eeirth’s Rotation and Intercompare the Techniques of observation and 
analysis (MERIT). The proposed techniques of observation included 
laser ranging and radio interferometry (Wilkins, 1980). As early as 1980 
(August-October) the MERIT Short Campaign was undertaken to test and 
develop the organizational arrangements that would be required for a 
realistic coordination and successful implementation of the MERIT Main 
Campaign which very successfully took place during the 14-month period 
of September 1, 1983, to October 31, 1984. 

The MERIT Working Group in collaboration with the Conventional 
Terrestrial Reference System (COTES) Working Group, decided to extend 
the objectives of the MERIT Main Campaign in order to include the 
preparation of a catalog with a precise and consistent set of station 
coordinates (Wilkins et al., 1986). Consistency of the station coordinates 
is achieved by accurately linking together the reference frames realized 
by each of the techniques involved. This is accomplished either 
through collocations or by estimating for each technique the diurnal 
differences of the earth rotation parameters. Towards achieving this 
goal, it was decided to have an intensive campaign, during which, in 
addition to other requirements, all the MERIT stations were asked to 
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observe as frequently as possible and in full capacity (ibid). The SLR 
technique, one of the techniques tested during the MERIT Main 
Campaign, reached its full potential. It would never have reached this 
potential if it weren’t for the NASA Crustal Dynamics Project whose SLR 
network after several years of buildup approached during the MERIT 
Main Campaign its full capacity (Coates et al., 1985). It was this peak in 
the operation of the SLR system that resulted in an extensive SLR 
simultaneous tracking which in turn sparked the initiation of the 
present study. Consequently, the SLR observations collected during the 
MERIT Main Campaign were used in this study and herein this data set 
is referred to as the MERIT Main Campaign (MMC) data set. 

According to the SLR organizational arrangements for project MERIT, 
each of the observing stations was obligated to submit its Full Rate (FR) 
observations to the Crustal Dynamics Information System (CDIS) located 
at Gioddard Space Flight Center (GSFC). The FR observations would 
have to be submitted within three months after their collection (Schutz, 
1983b). The MERIT FR data format is a character oriented format 
referred to as Seasat Decimal Format (SSD) (ibid.). 

The MMC data set is available to any investigator from the Crustal 
Dynamics Data Bank (CDDB). This data set was sent to us upon request 
in nine-track magnetic tapes. Each record in the tape is stored in SSD 
format and contains the observed range, the epoch of the observation 
and a number of indicators pertaining to the corrections that have been 
applied and to those that have yet to be applied (i.e., atmospheric 
refraction corrections, center of mass corrections, etc.). The observed 
ranges stored in those records are corrected for system delay, signal 
amplitude dependent effects and any other effects pertaining to the 
laser instrument (see Section 3.2). The center of mass correction has a 
negative sign, and therefore it should be applied to the computed range. 
Each record also contains the atmospheric refraction correction computed 
with the Marini and Murray formula together with the pressure, 
temperature and relative humidity recorded at the observing site. The 
meteorological data is included for the analysts who might prefer to 
compute the atmospheric refraction correction with a different formula 
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than that of Marini and Murray. The Marini and Murray formulation was 
suggested for use with the MMC data set as an extension of the MERIT 
standards (ibid.). 

At the time when this study was initiated the MMC data set was not 
available. Thusi at the initial stage of this work, we employed the 
Lageos laser ranges recorded during the last three months of the year 
1979 by stations 7114 (Owens Valley) and 7115 (Goldstone) located in 
California. The initial stage of this work was primarily devoted to the 
development and testing of the software necessary to edit the laser 
range observations and to generate the simultaneous range and 3RD 
observables which constitute the input to the geometric and 3RD 
methods respectively. Baseline estimation, however, was based solely on 
the MMC data set (see Chapter 4). 

The MMC data set is not a unique data set because for each month 
there have been several releases issued due to data problems or to 
missing data (Section 3.6.3). Table 3 shows the monthly releases of the 
MMC data set used in our study. These releases, besides having 

erroneous observations, also contain data records with unacceptable 
characters such as asterisks, plus and minus signs, etc. These records 
should not have been there since the received data set was supposed to 
have been preprocessed by the Bendix Field Engineering Corporation 
(Schutz, 1985). Editing this data set was a tedious process and a 
considerable amount of time was spent for this purpose. 
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Table 3. 

Monthly MERIT Releases 


Month 

Release 

Final 

Release 

Sep. 1983 

C 

yes 

Oct. 1983 

B 

no 

Nov. 1983 

E 

yes 

Dec* 1983 

E 

yes 

Jan. 1984 

E 

yes 

Feb. 1984 

C 

no 

Mar . 1984 

C 

no 

Apr. 1984 

D 

no 

May 1984 

D 

no 

Jun. 1984 

B 

no 

Jul. 1984 

B 

no 

Aug. 1984 

C 

no 

Sep. 1984 

B 

no 

Oct. 1984 

C 

yes 


3.5 DATA EDITING 

The brief description of the error sources affecting the SLR systems 
(see Section 3.2) reveals that erroneous timing of the returned laser 
pulse is possible, especially when the observations are made during 
daylight time with a single photon detection laser instrument (i.e., TLRS 
1, 2, 3 or 4) (Shawe et al., 1985). Erroneous timing is even worse for 
the laser instruments equipped with single stop Time Interval Units 
(TIU). With single stop TIU’s it is not possible to detect the returned 
laser pulse if a noise pulse with an energy level exceeding the stop 
discriminator threshold enters the receiver prior to the returned laser 
pulse. Multistop TIU's, on the other hand, have the potential to reduce 
erroneous detection substantially since they are designed to detect more 
than one returning pulse for each of the emitted pulses. The multistop 
TIU’s available on the market today have time resolution at the 
nanosecond level. Such low resolution makes them inadequate for use 
with centimeter-level accuracy laser instruments. Erroneous laser range 
observations may also result from operator errors, inadequate 

maintenance and from any other sources affecting the proper operation 
of the laser instrument. 
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The erroneous observations should be detected and rejected before 
the generation of the Simultaneous Range (SR) and the Simultaneous 
Range Difference (SRD) observables. In common practice, editing of the 
erroneous observations is incorporated in the final adjustment of the 
observations when station coordinates and baselines are estimated. In 
the present study, however, the editing of the laser ranges should 
precede the final adjustment simply because the SR and SRD observables 
are obtained through an interpolation of the laser ranges, and therefore, 
the presence of erroneous ranges will affect the entire set of the 
generated SR and SRD observables. Thus, effective generation of the 
SR and the SRD observables requires early editing of the observed laser 
ranges. 

3.5.1 Data Snooping Procedure 

Any kind of editing procedure requires a functional representation of 
the observed ranges. The estimation of the parameters involved in this 
representation allows the prediction of the observed ranges and 
therefore the estimation of the errors associated with each of those 

ranges. The difference of the estimated error vector (i.e., residual 

vector) from the true error vector accounts for the projection of the 
true error vector into the model space generated by the column space 
of A (see equation 3-12). This is the component of the true error that 

is lost in the estimation process, and therefore it cannot be recovered. 

Based on the statistical properties of the estimated component of the 
true error (i.e., residual), statistical tests may be derived to allow for 
the detection and rejection of erroneous observations. 

This section proceeds with the development of the statistical 
formulation necessary for the derivation and implementation of those 
tests. The next section contains a description and elaboration on how 
the base functions employed in the functional representation of the laser 
ranges should be chosen to effectively edit the laser range 
observations. 

Let’s consider the linear adjustment model employed in the 
functional representation of the observed laser ranges (Uotila, 1986) 
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(3-3) 


-e = AX - Lb 


where 

e = true error vector of the observations 
Lb = vector of the observed laser ranges 
X = true parameter vector 
A = design matrix of the experiment 

The minimum variance unbiased estimate of the parameter vector X has 
the following form (ibid.) 

X = (ATpA)-‘ATpLb (3-4) 

The true errors, under the null hypothesis, have a multivariate normal 
distribution with 0 mean and Ig variance-covariance matrix 

e ~ N(0,le) (3-5) 

where 

le “ = variance-covariance matrix of the observations (3-6) 

P = weight matrix of the observations 
<To = the a priori variemce of unit weight 

The unbiased estimate of the a priori variance of unit weight takes the 
following form (ibid.) 

^8 = j; (AX - Lb)Tp(AX - Lb) (3-7) 

where 

n - u = degrees of freedom of the adjustment 

The projection of the true error vector e into the orthogonal space 
complement to the model space constitutes the estimated error (i.e., 
residual) vector 

V = AX - Lb = [A(ATPA)-*ATP - I]Lb (3-8) 
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Substitution, in this equation, of the weight matrix P with the identity 
matrix leads to the following equation 

V = [A(ATA)-»AT - I] Lb (3-9) 

This substitution is well justified because the editing of the 
observations is carried out on a station-by-station basis and for each of 
the stations involved the variance-covariance matrix of their 
observations is assumed to have a diagonal form and equal diagonal 
elements. Thus, the weight matrix P, apart from a constant factor, is 
equal to the identity matrix. 

Since the true error vector e has zero mean, equation (3-3) takes 
the following form 

Lb = e + E(Lb) = e + AX (3-10) 

where the symbol E denotes the expected value. Substitution of 
equation (3-10) into (3-9) leads to the following equation 

V = Me (3-11) 

where 

M = [A(ATA)~^t - I] (3-12) 

The matrix M is a symmetric idempotent matrix with rank (n - u), and 
since MA = 0 it represents the orthogonal complement operator of the 
model space (i.e., projects any vector into the orthogonal complement of 
the model space) (Pope, 1976). 

With the help of equation (3—11) and the law of covariance 
propagation, the variance-covariance matrix of the residuals takes the 

following form 

Iv = M<rg (3-13) 

Iv = M9g (3-14) 
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The first factor of the matrix M (eq. 3-12) depends on the design 
matrix A or, in other words, on the chosen experiment. From an 

estimation point of view the experiment (i.e., design matrix) should be 
chosen to minimize the projection of the true observational error into 
the orthogonal complement of the model space (i.e., to minimize the effect 
of the true observational error on the predicted ranges). From an 

editing point of view in which we are interested, the M matrix should be 
nearly diagonal. This is preferred because each individual residual (v^ ) 
will be primarily affected by the corresponding true error e^ . This 
makes it easier to identify erroneous observations with one-dimensional 
residual testing (see next section). 

Using an independent set of analytic base functions for the 
functional representation of the laser ranges results in a severe 
limitation in regard to the choice of the model space. This limitation 
arises because any analytic base function can be approximated by a 
partial sum of monomials up to degree k. Thus, the space spanned by 
the monomials (1, t, . . . , t*^) closely resembles the model space spanned 
by any independent set of base functions. The basis, however, selected 
to span the model space will determine the conditioning of the normal 
matrix (i.e., A^A) and the distribution of the approximation errors in the 
interval of approximation. Choosing Chebychev polynomials as base 
functions results not only in a well-conditioned normal equations matrix 
but also in an even distribution of the residuals over the interval of 
approximation (see next section). 

Having a functional representation for the observed laser ranges we 
can proceed with the editing of those ranges using one-dimensional 
statistical testing of the residuals. This is accomplished with the 
one-dimensional data-snooping procedure originated by Baarda (1968). 
In this procedure, under the null hypothesis the true observational 
errors have a multivariate normal distribution with zero mean and Ig 
variance-covariance matrix (equation 3-5). Thus, under the null 
hypothesis Ho the residuals have a multivariate normal distribution 

V - N(0, Iv) (3-15) 
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with zero mean and ly variance-covariance matrix (equation (3-13)). 
This leads to a normal marginal distribution for each individual residual 
with zero mean and variance 

Vj ~ n(0,<r§^) (i = 1, . . . , n) (3-16) 

Thus> the original null hypothesis H© is now replaced by a sequence of 
null hypotheses Ho^ (i = 1, . . . , n): 

Ho : Vi ~ n(0, <r§i) (i = 1, . . . , n) (3-17) 

where 

Vi = the residual of the i^^' observation 

<TVi = <To '^577 (3-18) 

= a priori variance of unit weight 
iDi i = diagonal elements of the matrix M (equation 3-12) 

Under the sequence of the hypothesis H©. (i = 1, . . . , n) the statistic 
Wi = Vi/<x^^ has a standard normal distribution 

Wi ~ n(0, 1) (i = 1, . . . , n) 


The capital letter N in the above equations denotes multivariate normal 
distribution while the small letter n denotes one-dimensional normal 
distribution. 

The theoretical value (Tq in equation (3-18) is unknown and therefore 
the estimated 5© (i*e*> a posteriori variance of unit weight) is used to 
evaluate the W, statistic. This new (W^) statistic, under the sequence 
of null hypotheses H©. (i = 1, ... n), and as the degree of freedom 
becomes larger and larger tends to have a student’s t distribution 
(Pope, 1976) 
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(3-19) 


^ v^ 

W, = ^ 

where 

^Vi “ ®i 1 (3—20) 

t<o = Type I error = P (rejecting Ho I when Ho is true) 
n - u =: degrees of freedom 

a% =: a posteriori variance of unit weight (equation (3-7)) 

The critical region C of a statistical test based on the W< statistic takes 
the following form 

C = ( W^ : IwJ > t„-o, ) (3-21) 

while the Type I error becomes 



Therefore, the data-snooping procedure is carried out in the following 
steps: 

1. Choose a probability level (oto = 0.005) 

2. Take from a student t table the critical value (c) for rejection 
(e.g., c = ^n-u, oto, with n - u > 200 * 2.84) 

3. Compute the individual statistic 

4. Reject each observation Lt,. leading to > c 

If the observations are correlated a slightly different data snooping 
is required (Pope, 1976). The data-snooping procedure has been 
effectively applied in photogrammetry (Gruen, 1979), in deformation 
analysis (FIG Deformational Analysis Working Group, Heck 1982) and in 
the present study to edit laser range observations. 
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3.5.2 Effectiveness of the Data Snooping Procedure in Editing the 
Laser Range Observations 

Implementation of the data-snooping procedure requires a functional 
representation of the observed laser ranges. Before proceeding with 
the choice of the base functions, we should try to single out those 
properties of the base functions that would make the data-snooping 
procedure more effective. To accomplish this, we should realize that the 
data-snooping procedure is based on a sequence of one-dimensional 
tests carried out for each residual individually. These one-dimensional 
tests would be effective if the residuals are uncorrelated and evenly 
distributed over the interval of approximation. It is not possible to 
obtain uncorrelated residuals with the linear adjustment model (3-3); if 
this were possible it would imply that rank (M) = n, which is a 
contradiction since rank (M) = n - u (see equation (3-12)). Thus we can 
only look for base functions yielding residuals with reduced 
correlations. This, however, is not possible because a set of k 

independent analytic base functions can be uniquely mapped into the 
linear space spanned by the monomials up to degree k (i.e., 1, t, . . . , 
t**). Thus we are left only with the choice of base functions that yield 
an even distribution of the residuals over the interval of approximation 
and presumably a well-conditioned normal equations matrix. 

To understand how an uneven distribution of the residuals may 
result, we choose the monomials to represent the observed laser ranges 
in the interval [to, tg]. Since any such arbitrary interval can be 
transformed to the interval [-1, 1] by the change of variables. 



tg ^ t ^ tg 
-1 < T < 1 


(3-22) 


it is sufficient to examine the behavior of the monomials (1, t, t^, ... t**) 
in the interval [-1, 1]. In this interval each monomial assumes the same 
maximum absolute magnitude 1 at r = and the same absolute minimum 
magnitude 0 at r = 0. If the observed ranges would be approximated 
with the polynomial 
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R,(t) = St . t 


(3-23) 


where 

T = (1, T, T^, ... t‘‘)T 


(3-24) 


X = estimated monomial coefficient vector (equation (3-3)) 

A 

the errors in the parameter vector X will produce small residuals for 
small T (i.e., t near zero) and large residuals for t close to 1 or close 
to -1. Such uneven distribution will apparently place limits on the 
effectiveness of the data-snooping procedure. In addition to this, use 
of monomials as base functions gives rise to numerical problems 
associated with the inversion of a nearly-singular normal matrix when k 
is moderately large (Carnahan et al., 1969; Pavlis, 1982). These 
numerical problems further deteriorate the effectiveness of the 
data-snooping procedure. Thus, the monomials cannot be effectively 
used with the data-snooping procedure. 

To avoid uneven distribution of the residuals it seems reasonable to 
look for functions having evenly distributed extreme values of equal 
magnitude in the interval [-1, 1]. The Chebychev polynomials appear to 
be good candidates since their cos(nd) origin fulfills the above 
requirements. These polynomials are defined with the following 


equations 

T„(t) = cos(n0) ; n = 0, 1, ... k (3-25) 

where 

e = cos~Mt) (3-26) 

From this definition, one readily obtains 

To(t) = 1 (3-27) 

Ti(t) = T (3-28) 


and 
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T„(t) = 2t • T„-,(t) - T„-,(t) 


(3-29) 


The n real roots of the polynomial Tn(r) occur in the interval [-1, 1] at 
the points 

Ti = cos ; i = 1, ... k (3-30) 

Using equations (3-25) through (3-29) it is a matter of simple exercise 
to prove that the extreme values of Chebychev polynomials have the 
same absolute magnitude of (1) and are evenly distributed in the 
interval (-1| 1], In this interval the Chebychev polynomials are 

orthogonal with respect to a weighted integral operator with weight 
function w(t); 

w(t) = 1 /n/~ 1 - (3-31) 

We desire* however, that these functions be orthogonal with respect to 
summation as well. Fortunately, this is true but only if the summation 
is carried out over a specific set of points in the interval [-1, 1] 
(Pavlis, 1982) 


I Tiirt) Tj (Ti) 
i=0 


0 

m + 1 
2 

m + 1 


i * j 
i = j * 0 
i = J = 0 


(3-32) 


where rjt denote the roots of the polynomial T„+i(t) given by equation 
(3-30). Thus, in the formation of the normal equations matrix (A^A) 
there exists a tendency for cancellations among the products of 
different degree Chebychev polynomials according to the equation (3-32). 
This tendency prevents the numerical problems associated with the 
inversion of the normals and it Justifies the term "nearly orthogonal" 
often used when the Chebychev polynomials constitute the base 
functions in a least-squares adjustment. 
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The effectiveness of the data-snooping procedure to edit laser range 
observations is shown below for only two passes. Extensive 
experimentation, however, with sparse and dense data sets recorded by 
many different stations, indeed confirms that the results for these two 
passes are indicative for the overall performance of the data-snooping 
procedure. Pig. 1 shows the residuals, computed with equation (3-9), 
for two passes observed by stations 7114 (Owens Valley) and 7115 
(Goldstone) on October 31, 1979. These residuals indicate that erroneous 
observations with blunders as large as 350 m do exist in the original 
data set. Rejection of the erroneous observations is carried out by the 
data-snooping procedure (see Fig. 2). This figure shows the 
distribution and the magnitude of the residuals after the application of 
the data-snooping procedure. It is evident that a rejection of less than 
10% of the the observations not only eliminates the blunders and makes 
random the residuals but also reduces the RMS from about 46 m down to 
0.11 m. 

A close inspection of Pig. 1 and 2 also reveals that observations 
having 20 m residuals not only survived in the data-snooping process 
but also reduced their residuals down to the 0.10 m level. This not 
only demonstrates the dependence of the residuals on the number and 
the magnitude of the blunders affecting the observations but also shows 
how difficult it would be to detect and reject erroneous observations by 
testing the individual residuals v^ alone. The data-snooping process 
overcomes this difficulty by testing the normalized residual (i.e., 

instead. Furthermore, Fig. 2 shows that about three percent of 
the residuals have an absolute magnitude greater than 0.40 m while the 
rest of them are randomly distributed about zero with an RMS of 0.12 m. 
Thus, it was decided to edit out the observations with residuals of 
absolute magnitude greater than 0.40 m. A plot of the residuals for the 
remaining observations is shown in Fig. 3. These residuals have a 
random distribution about zero with an RMS ranging from 0.10 m to 0.12 
m. These RMS values are consistent with the expected accuracy of the 
observing stations thereby confirming our claim that these observations 
form a clean (i.e., no presence of blunders) data set. Although the 
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Fig. 1 Residuals of Chebychev interpolation (before data snooping). 
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Fig. 1 (cont’d) 
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Fig. 2 Residuals of Chebychev interpolation (after data snooping). 
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Pig. 2 (cont’d) 


data-snooping procedure is very effective, it is relatively slow and 
therefore very expensive since five to six iterations are needed for the 
editing process to finish. Alternative ways to speed up this process 
have been developed by various investigators (Pope, 1976; Gruen, 1979). 
Finally, we should keep in mind that constant and time dependent 
biases, no matter how large, cannot be detected by the described 
procedure simply because if such biases do exist they will be absorbed 
by the recovered coefficients of the Chebychev polynomials. 


3.6 GENERATION OP SIMULTANEOUS RANGES AND SIMULTANEOUS RANGE 

DIFFERENCES 

Since Lageos is a passive satellite, it is quite unlikely for the 
coobserving stations to record strictly simultaneous observations even if 
the same part of the the Lageos orbit is coobserved. This happens 
because for each observing site the tracking starts at a different epoch, 
each laser instrument has a different repetition rate and last but not 
least there will always be synchronization errors among the coobserving 
sites. Implementation of the geometric and the SRD methods requires 
strict simultaneity, and therefore, an interpolation of the observed laser 
ranges is necessary. 

Simultaneous observations for the geometric solution are obtained by 
first identifying passes continuously coobserved (i.e., data gaps smaller 
than 60 seconds) by four or more stations (see Sections 3.6.2 and 3.6.3). 
For each of those passes the station with the least number of 
observations is identified. At its observing epochs simultaneous 
observations for all the remaining stations are generated through an 
interpolation. 

Simultaneous Range-Differences for the SRD method are obtained by 
dividing the observing stations into pairs with quasi-simultaneous 
observations. For each of these pairs the station with the least number 
of observations is identified. Subsequently, at its observing epochs 
interpolated ranges for the alternate station are generated. The SRD 
observables are finally obtained by subtracting the actually observed 
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Fig. 3 (coni'd) 


ranges of the station with the least number of observations from the 
corresponding interpolated ranges of the alternate station. 

ThereforOi it is essential for the success of this study to select an 
interpolation method that is capable of generating laser ranges with an 
accuracy compatible to that of the observations. 

3.6.1 Chebychev Polynomials and Spline Functions in the Context of 
Global and Piecewise Interpolation 

A survey of the interpolation methods shows that these methods may be 
divided into two basic categories: 

- the global interpolation methods, and 

- the piecewise interpolation methods. 

With the global interpolation methods, a function is approximated over 
the entire interval of approximation by the same linear combination of a 
selected set of base functions. The function being approximated may be 
known either analytically or quantitatively at a small number of base 
points. The latter case depicts the situation in the present study since 
the ranges to the satellite, apart from measurement errors, are known 
only at each of its observing epochs. With a piecewise interpolation 
method, a specific function whose values are given at a specified set of 
base i>oints is approximated by dividing the base points into successive 
subsets, each of which contains two, three or more base points. Within 
each subset the function is approximated by a different linear 
combination of (possibly) different base functions. Boundary conditions 
are imposed on the common points of adjacent subsets to make the 
interpolating function continuous with (possibly) continuous first- 
and/or second-order derivatives over the interval of approximation. 
Thus, with a piecewise interpolation method one obtains a continuous 
interpolating function consisting of pieces, each of which is composed 
from a different linear combination of (possibly) different base 
functions. 

The approximating function of a global interpolation method, on the 
other hand, generates approximate values with a relatively strong 
dependence on the values the approximated function assumes at each 
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base point* This is a desired property for our study because we would 
like to obtain interpolated ranges that reflect to the highest degree of 
accuracy the overall information inherent in the actually observed 
ranges. Furthermore, the effect of the gaps (i.e., distances between 
successive base points) should be effectively controlled, and if possible, 
kept down below the noise level of the observations. In global 
interpolation methods these effects of unevenly distributed gaps are 
uniformly distributed over the interval of approximation provided that 
the gaps are not large enough to corrupt the effectiveness of the global 
interi>olating function (see Section 3.6.2). With global interpolators the 
generated approximate values exhibit strong correlations because the 
same linear combination of base functions is used to generate 
approximate values over the entire interval of approximation. Moreover, 
the closer the approximated values the stronger the correlations are. 
Strong correlations are not welcomed in the present study since the 
generated SR and SRD observables will be considered uncorrelated in 
the final adjustment when the station coordinates and baselines will be 
estimated (see Sections 3.6.2, 4.3 and 4.4). The effectiveness of the 

global interpolation methods is largely dependent not only on the choice 
of the base functions but also on their implementation to "best" 
represent the given set of data points. 

The base functions most often encountered in practice are the 
monomials, the Chebychev polynomials, the Fourier series, the 
exponentials, etc. (Carnahan et al., 1969; Davis, 1975). Linear 
combinations of either monomials or Chebychev polynomials are by far 
the most important and most popular approximating functions. These 
polynomials are easily operated on by addition, multiplication, 
integration, differentiation, scaling and shifting and more importantly 
they are closed with respect to any of these operations. Other base 
functions also possess some or all of the above properties, and therefore 
the polynomials would not be so important if it were not for the 
Weierstrass approximation and uniform approximation theorems (Davis, 
1975; pp. 24 and 107). 
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The choice of the base functions depends on the behavior of the 
function being approximated. For instance, functions with certain 
periodicities can be best approximated with Fourier base functions 
sin(kT), cos(kr), k = 1, ... n, while functions with exponential behavior 
are best described with exponential base functions. A combination of 
Fourier series and exponentials would also be effectively used if the 
behavior of the approximated function exhibits such a pattern. In the 
present study, since the laser range observations are affected by 
periodic and secular perturbations caused mainly by the gravity field it 
is only fair to choose for their representation periodic base functions 
supplemented with monomials of degree zero, one, and possibly two. 
Chebychev polynomials exhibit such a behavior not only because of their 
cos(n0) origin (see Section 3.5.2) but also because the Chebychev 
polynomials of degree zero and one coincide with the monomials of the 
same degrees. In addition to this, the optimal properties of the 

Chebychev polynomials (see Section 3.5.2) make them ideal base functions 
for an effective global representation of the observed laser ranges. 

Having chosen a set of base functions, their linear combination 
should be determined to "best" approximate the function implied by the 
given set of data (i.e., function of the observed laser ranges). The best 
representation of this function is determined on the basis of a chosen 
criterion. Such a criterion may be chosen to either reproduce the 
function at its base points or to reproduce the function and its 
derivatives at a given point. We may also choose to either minimize the 
maximum error of the approximation (i.e., minimax principle) or to 
minimize a weighted sum of the squares of the residuals at the base 
points (i.e., least squares approximation). 

Reproducing the function at n base points results in an 
interpolating polynomial of degree (n-1). Thus approximation of the 
laser ranges with an interpolating polynomial would result in a 
polynomial of degree ranging anywhere between 200 to 12000. Use of 
such a high-degree polynomial is completely out of the question. 
Reproducing the function and its derivatives at a given point or 
minimizing the maximum error of approximation are also ruled out 
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because for the former criterion the derivatives needed are not available 
while for the latter the base points should coincide with the roots of 
the minimax polynomial. Such a choice is not feasible in our study 
simply because we have no control over when the laser ranges will be 
recorded. However | if such a choice would have been possible it would 
result in a minimax polynomial of very high degree} and therefore it 
would have been ruled out again. Minimizing a weighted sum of the 
squares of the residuals (i.e., least squares principle) constitutes a very 
good alternative not only because we can statistically select a relatively 
low degree for the approximating polynomial but also because the least 
squares do not reproduce the observations. The former property 

prevents instability problems usually associated with the polynomial 
interpolation} while the latter is desired because the observed laser 

ranges are always contaminated by measurement errors. Therefore, the 
functional representation of the observed laser ranges with Chebychev 
polynomials whose coefficients are estimated with a least squares 
adjustment} constitutes an alternative having many of the desired 
properties necessary for a successful interpolation (see Section 3.6.2). 

The piecewiae interpolation methods generate approximate values that 
are sensitive to the values of the approximated function in the 
neighborhood of the interpolating point but largely insensitive to the 
values of this function a little farther away from that point. 

Furthermore, the effects of the gaps in the piecewise interpolation 
methods are not uniformly distributed over the interval of 

approximation. Thus, approximation errors committed closer to big gaps 
are substantially larger than those committed closer to smaller gaps (see 
Section 3.6.2). Such behavior of the approximation errors cannot be 
easily controlled because their functional dependence on the magnitude 
and the distribution of the gaps is not known. Direct evaluation of this 
functional dependence requires laser ranges between adjacent 
observations which of course are not available. The inability to 

effectively control the effect of the gaps in the piecewise interpolation 
methods constitutes a major drawback when these methods are compared 
to the global interpolation methods (see Section 3.6.2). 
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The generated approximate values via piecewise interpolation 
methods are largely uncorrelated because a different linear combination 
of a (possibly different) set of base functions is used to generate the 
approximate values between different pairs of successive observations. 
This property differentiates the piecewise interpolation methods from the 
global interpolation methods and constitutes a desired property for the 
present study because the generated SR and SRD observables will be 
considered uncorrelated in the final least squares adjustment when 
station coordinates and baselines are estimated (see Sections 4.3 and 
4.4). 

With the piecewise interpolation methods a relatively small number of 
base functions (usually three to five) is needed to represent the data 
between adjacent data points. The choice of the base function should 
be made along the same lines discussed in the global interpolation 
methods. The coefficients of the different linear combinations of the 
(possibly) different set of base functions for each subset of adjacent 
base points are determined by imposing boundary conditions to 
reproduce the functional values (i.e., observed ranges) and to obtain a 
continuous interpolating function with continuous first and/or second 
derivatives. One may also choose not to reproduce the observed ranges 
but rather to introduce a weight function that would reflect a desired 
relation between the predicted and the actually observed laser ranges. 
The introduced weight function may be derived according to the 
measurement errors and to the distribution of the gaps in the 
neighborhood of the base points (i.e., observing epochs). Use of weight 
functions with cubic splines leads to the weighted cubic splines 
interpolation (Spath, 1974). Estimating the weight function is a difficult 
task, and this function may not be valid for different cases. Therefore 
it was considered appropriate not to introduce any weight function but 
rather to reproduce the functional values (i.e., observed ranges) at each 
base point (i.e., observing epoch). This would enable the evaluation of 
the overall performance of the piecewise methods in regard to their 
ability to generate good approximate values when relatively large gaps 
exhibit an uneven distribution over the interval of approximation (see 
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Section 3.6.2). For this evaluation we made use of the easy-to-handle 
cubic splines. 

A cubic spline s(t) where t is in the interval [tj, t„] is defined as a 
set of third-degree polynomials, each defined in the interval [tj,, t|<+i ] 
(k = 1, ... n-1). These polynomials are joined at the base points (t|() so 
that the resulting cubic spline is twice differentiable at each base point. 
A cubic spline defined on an interval with n base points consists of 
(n-1) third-degree polynomials. Thus, its unique determination requires 
4(n-l) independent sets of equations. The requirement to have a twice 
differentiable cubic spline in the interval [ti, t„] introduces a set of 
2(n-2) independent equations which results from the continuity 
conditions required for the existence of the first and second derivatives 
at the base points tj through t„_,. The requirement to reproduce the 
functional values (i.e., observed laser ranges) at the base points (i.e., 
observing epochs) introduces another set of [2(n-2) + 2] independent 
equations bringing up to (4n-6) the total number of independent 
equations. The two additional equations, necessary for the unique 
determination of a cubic spline, are obtained from the boundary 
conditions specified for both ends of the approximation interval (i.e., at 
the base points tj and t„). Since the observed laser ranges vary 
slowly within a few seconds of time, we have adopted in the present 
study the following conditions: 

s"(tx) = s"(t 2 ) (3-33) 

s''(t„) = s"(tn-i) (3-34) 

Such an arbitrary choice influences the results very slightly in the 
neighborhood of the end points (Spath, 1974; Pavlis, 1982). The 
independent set of equations used to determine the (4n-4) coefficients of 
a cubic spline is given in (Spath, 1974; Pavlis, 1982). The next section 
contains an evaluation of the relative performance of the least squares 
Chebychev and cubic spline interpolators. This evaluation is based on 
their ability to effectively interpolate the observed laser ranges. 
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3.6.2 Chebychev Polynomials vs. Cubic Spline Functions in the 
Functional Representation of Laser Ranges 
Least-squares approximation of the laser ranges with Chebychev 
polynomials requires a priori knowledge for the degree of the resulting 
algebraic polynomial 

Ra(T) = o<o + o{iTi(r) + . . . + o(|<T|((t) (3—35) 

The degree of the algebraic polynomial Ra(r) coincides with the degree 
of the Chebychev base function T|{(t). This degree should be 

chosen to represent the data "sufficiently" in the sense that an 
extension of equation (3-35) to include higher-degree Chebychev 
polynomials will represent the data with the same accuracy. There is a 
limit as to what degree we can go to because after a certain degree has 
been reached instability problems will deteriorate the solution 
substantially. With Chebychev polynomials there exists a wide range of 
degrees that can be used to represent the data with the same accuracy. 
This is possible because Chebychev polynomials are not seriously 
affected by instability problems (see Section 3.5.2). This property will 
be very useful in our study because in the generation of the SR and 
SRD observables it reduces the computing time substantially (see Section 
3.6.3). 

The lowest-degree Chebychev polynomial that sufficiently represents 
the available data is conveniently determined through statistical testing. 
To construct this test it was assumed that the coefficients of the 
equations (3-35) were estimated together with the weighted sum of the 
squares of the residuals 

Rk = (VTpV)k (3-36) 

through a least-squares adjustment. At this point we would like to 
statistically test if the coefficient (o(|() is significantly different from 
zero. To accomplish this, we perform another adjustment supplemented 
with the constraint that the coefficient (ofk) is zero. If in this 
adjustment the computed weighted sum of squares of the residuals is 
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Ri(_i, then the weighted sum 


Rc - Rk-i 


- Rk 


(3-37) 


caused by the constraint (o<,, : 0) is distributed as <r*xi I = <y*N^(0, 1)] 
independently from R|< (Hamilton, 1964). Thus, the percentage change of 
the weighted sum of the residuals 


Rr Rk— 1 Rfc 

P = -(n-k-1) = -(n-k-1) 


Xi 


Xn-k-i 


“ Rl, n-k-1 - 


(3-38) 

caused by the constraint (of^ = 0) is distributed as t*_k_i. Having P 
and its distribution it is possible to set up a hypothesis test based on 
the percentage change of the weighted sum of the residuals 


Hq: = 0 ( = : not significant change in P) 

Hi: o(|, * 0 ( = : significant change in P) 

The critical region of the test is obtained through the following equation 

P(Fi, n-k-1 > Per) = o( P(|t„-k-il > (Per)’^) = | (3"39) 

where <x is the significance level (i.e.. Type I error) of the test. 
Specifying a value for o< of 0.01 (1%) we can determine P^r for n > 200 
with the help of a t table and equation (3-39) (DeGroot, 1975) 

(Per)^ = 2.576 

resulting in the critical region 

C = {P £ R: P > (2.576)*} (3-40) 

Thus the hypothesis Ho is rejected if 

P > (2.576)* (3-41) 

thereby suggesting that a higher degree is needed at the 1% 
significance level. This process continues up to the degree for which 
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Ho cannot be rejected. With this test at our disposal we have the 
means to determine the least squares Chebychev polynomial that 
^'sufficiently” represents the observed laser ranges. Therefore, we can 
proceed with the relative evaluation of the cubic spline and Chebychev 
interpolators. This evaluation is based on the ability of these 
interpolators to "adequately” represent the laser ranges. 

The relative performance of these interpolators is based on a 
quantitative analysis of the orbit residuals obtained with two identical 
semidynamic orbit adjustments. The SRD observables input to those 
adjustment were generated through the cubic spline and the Chebychev 
interpolators using the edited observations of the station pair 7114 and 
7115, These observations correspond to the same two passes whose 
Chebychev residuals are shown in Fig. 3. The generated cubic spline 
and Chebychev SRD^s were processed with the help of 

- Lageos initial state vectors predicted, through numerical 
integration, for the starting epochs of both passes 

- the shortened version of the DE/LE200 lunar planetary ephemeris 
file covering the time span of the observations, and 

- the coordinates of the pole with the variations in UTl at each 
observing epoch 

through two separate orbit least squares adjustments implemented by 
the GEOSPP software. This software was developed by Pavlis (1982) and 
during the course of this study was not only modified to comply with 
the MERIT standards but also was corrected to successfully operate in 
the real data environment (see Chapter 2), The TRF and the CRF frames 
for each of the above adjustments are realized through the implicit 
constraints discussed in Chapter 2 supplemented with the following 
weighted constraints (see Section 4,1): 

- coordinates of station 7114 were held fixed (i.e,, ax = rry - az 
0,0001 m) 

“ coordinates of station 7115 were moderately weighted (i,e,, ax = 
ay = az = 10.0 m) 

- initial state vectors were moderately weighted (i,e,, ax = ay - az 
= 20 m; ax = ay = az = 0,02m/s) 
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The relative evaluation of the cubic spline and Chebychev interpolators 
is based on the comparison of the orbit residuals obtained with the 
cubic spline and Chebychev SRD observables. The cubic spline and 
Chebychev orbit residuals for the same two passes depicted in Fig. 3 
are shown in Figs. 4 and 6. The generation of the cubic spline SRD’s 
and the Chebychev SRD’s is based on the interpolation of the station 
having the larger number of observations. 

A comparison of Figs. 4 and 5 shows a strong correlation between 
the cubic spline orbit residuals and the data gaps of station 7115. 
However, this strong correlation exists only for the pass whose 
successive data gaps reach a maximum value of 60 seconds (see left 
plots of Figs. 4 and 5). For the other pass, whose gaps are not larger 
than 10 seconds, the correlation of the cubic spline orbit residuals with 
the successive data gaps is not strongly pronounced. This happens 
because the small magnitude data gaps are evenly distributed over the 
entire pass (see right plots of Figs. 4 and 5). On the contrary, the 
correlation of the Chebychev orbit residuals with the successive data 
gaps of station 7115 is very low (compare Figs. 6 and 5). Furthermore, 
it is quite clear that the cubic spline residuals are noisier than their 
Chebychev counterparts. In fact, the RMS of the cubic spline residuals 
for the passes shown in Figs. 4 and 6 is 0.41 m and 0.23 m respectively, 
while for the same two passes the RMS of the Chebychev orbit residuals 
is 0.12 m and 0.11 m respectively. The noisy behavior and the high 
RMS of the cubic spline orbit residuals is traced to the piecewise nature 
of the cubic spline interpolation. This assessment is confirmed by 
computing the differences between the cubic spline and the Chebychev 
SRD’s. These differences are shown in Fig. 7, and they seem to exhibit 
an almost identical behavior with the cubic spline orbit residuals (see 
Fig. 4), thereby confirming the anticipated fact, namely, that these 
residuals are primarily caused by the cubic spline interpolator. This 
result is not surprising since the piecewise nature of the cubic spline 
interpolator makes it sensitive to the magnitude and distribution of the 
data gaps. On the contrary, the global nature of the Chebychev 
interpolator makes it insensitive to the magnitude and the distribution 
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Fig. 4 Orbit residuals (cubic spline SRD’s). 



Fig# 5 Successive data gaps for station 7115# 








of the data gaps provided that these gaps are not larger than 60 
seconds. On the other hand, the 0.11 m and 0.12 m RMS of the 
Chebychev orbit residuals is consistent with the expected accuracy of 
stations 7114 and 7115 respectively. Thus one can safely assume that 
the errors caused in the generation of the Chebychev SRD^s by the 
up-to-60-second data gaps are not larger than the noise level of the 
observations. Extensive experimentation with dense and sparse data 
sets has revealed that gaps larger than 60 seconds tend to be several 
minutes long, thereby implying that the corresponding station ceased to 
observe due to calibration, weather problems, etc. (see Section 3.6.3). 
Thus, from now on we say that a station observes continuously if and 
only if successive data gaps are not larger than 60 seconds. 
Furthermore, interpolation is performed only over time intervals with 
"continuous^' coverage of observations. 

The above discussion demonstrates that interpolation of the laser 
ranges with Chebychev polynomials is superior to that of the cubic 
splines. Thus, for the MMC data set all the SR and SRD observables 
have been generated through a Chebychev least squares approximation. 
The interpolation is performed only over observing intervals with gaps 
smaller than 60 seconds (see next section). 


3.6.3 Data Selection for the Gteneration of the Simultaneous Ranges and 
the Simultaneous Range Differences 

In the present study a considerable amount of time and effort was 
devoted to generate the Simultaneous Range (SR) and the Simultaneous 
Range Difference (SRD) observables. These two observables constitute 
the input to the geometric and to the SRD methods respectively. A 
geometric solution, apart from special cases, is possible if at least six 
satellite positions are being coobserved by at least four stations 
generally distributed in space (i.e., not lying on the same plane). If the 
ground stations either form a plane (or close to forming a plane), a 
geometric solution is possible only if six or more stations are involved 
and each satellite position is coobserved by at least four stations (see 
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Section 2,1.3). An SRD solution, apart from special cases, is possible if 
at least one pair of stations is coobserving (see Section 2.2.6). 

The only data set having the potential to best fulfill the above 
reqirements is the MMC SLR data set (see Section 3.4). This data set 
contains the satellite laser range observations collected during the 
MERIT Main Campaign. Presently, a large amount of simultaneous laser 
range observations are recorded in the WEGENER/MEDLAS project. 

The MMC data set was received upon request from the CDDR data 
bank on nine-track magnetic tapes, each containing a month’s worth of 
observations collected either by all or by some of the MERIT stations. 
In each tape the information relevant to a specific observation is stored 
in records in the Seasat Decimal (SSD) format (see Section 3.4). Each 
record contains the observed laser range, the epoch of the observation, 
systematic corrections and indicators for the corrections that have been 
applied and for those that have yet to be applied (Schutz, 1983b). Some 
observations, however, were missing from the received tapes due mainly 
to data problems and delays. Thus, it was necessary at first to identify 
for each of the received tapes the ID’s of the observing stations, the 
number of passes per station and the number of observations per pass. 
This information was compared against the same information published in 
the monthly SLR reports issued for the entire MERIT Main Campaign by 
the Center of Space Research (CSR) at the University of Texas (UTX). 
These monthly reports contained the number of passes together with the 
total number of observations recorded by each of the stations involved. 
Any observations missing from the received tapes were obtained through 
a subsequent request issued only if the utilization of those observations 
was considered critical for the successful implementation of either the 
geometric or the SRD methods. Upon this request, new tapes identified 
with the release letters A, B, C, etc. were received from the CDDB data 
bank. These tapes contained the missing observations together with 
those received in previous releases. 

When the software necessary to process the laser range 
observations was ready, it was decided not to wait any longer but 
rather to use whatever releases were available up to that date (i.e., 
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Spring 1985), These releases are listed in Table 3, and they were 
finally employed in our investigation. For each of those releases the 
passes containing arcs coobserved by two or more stations were 
identified and isolated. In the context of the present study a pass is 
defined as a satellite passage whose blind spots are shorter than 2.2 
hours. With blind spots we designate the periods during which the 
satellite was not observed by any of the stations involved. An arc, on 
the other hand, denotes an observing period during which the 
successive data gaps are shorter than 60 seconds. An arc overlap of a 
Lageos pass recorded by American stations is shown in Figure 8 (left 
plot). 

The abscissae in these two plots denote the observing epochs 
relative to the starting epoch of the pass shown on the top of these 
plots. The starting epoch (84082442132) shown on the top of the h'f! 
plot reads 84 (1984), 08 (August), 24 (-th day), 4 (hours), 21 (minutes) 
and 32 (seconds). The ordinates on the right-hand side designate tri?' 
numeric station ID's (see Table 2 in Section 3.2.4). The ordinates on the 
left-hand side of these plots designate in ascending order the observing 
sequence of the stations involved. The solid horizontal lines indicate 
that for the period they cover the stations whose numerical ID is shown 

on the right-hand side of this line have data gaps shorter than 60 

seconds. The dotted horizontal lines indicate that for the period the\^ 
cover the corresponding station was not observing. The number above 
the starting point of each solid line denotes the number of observations 
whose successive data gaps are shorter than 60 seconds. Fig. 8 

confirms the assertion made in the previous section, namely, if the gaps 
are larger than 60 seconds they tend to be several minutes long. It is 
evident from the left plot of Fig. 8 that all of the stations invoK'ed 

except 7086 have a good continuous coverage, thereby making it possible 
to effectively interpolate their ranges. There is no need to interpolate 
the ranges of station 7086 since this station has recorded the least 
number of observations (see Section 3.6). The observing pattern shown 
in the left plot of Fig. 8 is indicative for the performance of the 
American stations during the MERIT Main Campaign. On the contrary. 
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Fig. 8 Arc overlap of two Lageos paBsea recorded by American stations (left plot) and 
European stations (right plot) respectively. 



the performance of the European stations is not quite as good. This is 
seen in the right plot of Fig. 8. This plot confirms the extended 
overlap expected among the visibility regions of the European stations. 
In spite of this, we can hardly interpolate the laser ranges for any of 
the coobserving stations simply because each of the arcs involved 
contains a small number of observations. Thus the question as to how 
many continuous observations (i.e., observations with successive data 
gaps shorter than 60 seconds) are enough for an effective interpolation 
is examined next. This question will be investigated by analyzing the 
errors committed in the recovery of ground truth observations. The 
recovery errors are computed by taking out a subset of observations, 
referred to as ground truth observations, in such a way that the data 
gaps of the observations left remain shorter than 60 seconds. The 
observations left are then used to deteiaiiine a least-squares Chobychev 
polynomial which is subsequently employed to recover the ground truth 
observations, that is, the observations not considered in its 
determination. Finally, by subtracting the recovered ground truth 
observations from the actually observed ones we obtain the re'cnvf^ry 
errors mentioned above. A sample plot of such errors for station 7210 
is shown in Fig. 9 together with the distribution of the ground truth 
points used in the computation of those errors (i.e., left- and 
right-hand-side plots respectively). 

These errors were computed for 285 ground truth observations 
obtained from a total of 1427 available observations by taking out every 
fourth observation. The ones left (i.e., 1142 all together) formed the 
basis to determine the least-squares Chebychev polynomial which was in 
turn employed to recover the ground truth observations. The abscissae 
for both of the plots shown in Fig. 9 designate the epochs of the 
ground truth observations relative to 0*^ UT of the day shown on Lhe 
top of the figure. The ordinates of the left plot designate the recovery 
errors (i.e., Range - Hint.) while the ordinates of the right plot 
designate the laser ranges after they have been scaled and shifted for 
plotting purposes. The RMS of these errors is 0.03 m which is exactly 
equal to the expected precision of station 7210. 
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Fig. 9 Recovery errors and distribution of ground truth observations (dense data). 






Using a different arc observed by the same station and in the same 
month we a^ain have computed the recovery errors for -17 ground truth 
observations obtained from a total of 283 observations by taking out 
every sixth observation. These errors, together with their distribiitions, 
are shown in Fig. 10, 

The RMS of these errors is about 23% worse than the expected 
precision of station 7210 (Analysis of Lageos Range Data, August 1985), 
Extensive experimentation with American stations has indeed confirmed 
that on the basis of less than 500 observations the interpolated ranges 
are affected by errors which are 15% to 20% worse than the noise; level 
of the observations. On the contrary, using more than 500 observations 
the errors caused by the interpolation hardly ever reach the noise level 
of the observations. 

The distribution of the ground truth points shown in Figs. 9 and 10 
reflects for each arc the distribution of the available observations 
because the ground truth observations were obtained from the available 
observations by taking out either every sixth or every foiirth 
observation. With this in mind a closer inspection of Figs, 9 and 10 

cleaidy reveals that small recovery errors are distributed around the 
denser parts of these two arcs while for the denser arc shown In P:g. 9 
the recovery errors are considerably smaller. Although the small 
recovery errors are distributed around the denser parts of thesf' two 
arcs, there exist in those parts recovery errors having an absolute 
magnitude of 0.10 in which is about three times larger than the eAX)ec:led 
accuracy of station 7210, Thus, the question arises as to w^hether these 
large recovery errors are caused by the procedure used to compute 
them or by errors affecting the actually observed ranges. 

To examine this we have plotted in Fig. 11 the Chebychev residuals 
for both arcs shown in Figs, 9 and 10. These Chebychev residuals have 
been computed by using all of the 1142 and 283 available observations 
for both arcs. Comparison of the recovery errors shown in Figs, 9 and 
10 with the corresponding Chebychev residuals shown in Fig. 11 clearly 
shows their high correlation and the equality of their RMS values as 
well. Moreover, th^; large recovery errors seem to have their 
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Fig. 10 Recovery errors and distribution of ground truth observations (sparse data) 








counterparts in the Chebychev residuals thereby suggesting that the 
large recovery errors have originated from errors affecting the actually 
observed ranges. More importantly the equality of the RMS values 

between the recovery errors and the Chebychev residuals suggests that 
analysis of the Chebychev residuals will give a good qualitative measure 
for the accuracy of the interpolated ranges. This important result used 
in the present study to control the quality of the interpolated ranges 
has also been confirmed from the analysis of the orbit residuals (see 
Section 4.4). Thus, we can safely state that interpolated laser ranges 
based on a relatively dense data set (i.e., with more than 500 
observations) will be affected on the average by approximation errors 
that are smaller than the noise level of the observation. Since the 
Chebychev polynomials are nearly orthogonal, deterioration due to 
numerical instability ceases to exist as the number of observations 
increases 11,000, 12,000 or more. Thus, since the recovery errors are 
smaller for denser data sets it is preferable to interpolate dense rather 
than sparse data sets provided that at least 500 observations with gaps 
shorter than 60 seconds are available. 

Inspection of the left and right plots of Fig. 8 reveals that the 
American stations fulfill the above requirements (i.e., to have more than 
500 observations with gaps shorter than 60 seconds) while the European 
stations hardly ever fulfill such a requirement. The American station 
7086, on the other hand, seems to have experienced several problems 
during the pass shown in the left plot of Fig. 8. The observing 
pattern, however, of this station was more or less the same for the 
entire MERIT campaign thereby making its observation inappropriate for 
interpolation. Fortunately, out of sdl the American stations used in this 
investigation only two of them, station 7086 (Texas) and 7112 (Colorado), 
seem to consistently have so many interruptions within all of their 
observed passes. More specifically, the interruptions of station 7086 are 
traced to the need of calibrating the laser within any of its observed 
passes. Furthermore, for most of the passes coobserved by four or 
more American stations, there exist at least three stations with dense 
enough observations to be effectively interpolated. This is the 
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necessary requirement only for the geometric solution since the required 
four-station events are obtained by interpolating the observations from 
three stations at the observing epochs of the fourth station for which 
the actual ranges are used. For an SRD solution, the observations of 
only one station are interpolated at the observing epochs of the 
alternate station. 

The European stations, on the other hand, seem to have experienced 
several interruptions for most of their observed passes. These 
interruptions have forced most of them to record relatively sparse 
observations. Therefore, it is difficult to find passes having arcs 
coobserved by four or more stations out of which three have dense 
enough observations that can be effectively interpolated. For instance, 
for the arc overlap shown in the right plot of Fig. 8 the observations 
from stations 7839, 7834 and possibly from 7840 can be effectively 
interpolated while the observations from any of the remaining stations 
are inappropriate for interpolation. Since this is the case for most of 
the passes coobserved by the European stations, it was considered 
appropriate to drop these stations from their implementation in the 
geometric solution. However, simultaneous observations collected in 
(WEGENER/MEDLAS) project have the potential to effectively implement 
the SRD method. 

In other parts of the world there may exist two or three stations 
coobserving. Overlap, however, for four or more stations is quite 

unlikely to occur anywhere else except over North America and Europe 
because only in these two parts of the world there exists a large 
number of operating laser ranging stations. This and Lageos* high 
altitude orbit are the main reasons why most of the passes recorded 
during the MERIT Main Campaign by American stations were coobseved 
by at least two of them. This resulted in a large number of American 
station pairs with quasi-simultaneous observations, and therefore it was 
possible to generate a large enough number of SRD observables that 
have led to a steady state response of the semidynamic solutions (see 
Chapter 4). Since the purpose of the present study is not to compute 
all possible baselines but rather to study the performance of the SRD 
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and geometric methods in baseline determinations, it was considered 
appropriate to employ the laser ranges recorded only by the American 
stations for both the SRD and the geometric methods. Fig. 12 shows the 
locations for those stations. 

Since for each station the observing pattern in terms of data gaps 
and observational density is relatively homogeneous from pass to pass, 
weVe shown in Fig. 13 three additional arc overlaps involving stations 
that are shown in Fig. 12 but not in Fig. 8. A comparison of Figs. 8 
and 13 confirms our assertion that the observing patterns of stations 
7109, 7110 and 7122 which are involved in more than one pass are 
homogeneous from pass to pass. Furthermore, the RMS of the 
Chebychev residuals for each of the stations involved shows very small 
fluctuations from arc to arc. This, however, is not true for stations 
such as 7122 that were upgraded during the MERIT Main Campaign. 
Table 4 lists the monthly mean RMS values of the Chebychev residuals 
obtained by interpolating the edited with the data-snooping-procedure 
laser range observations. The RMS values listed in Table 4 are in close 
agreement with those computed by the CSR at UTX listed in Table 2 in 
Section 3.2. This close agreement confirms once again the effectiveness 
of the data-Snooping procedure to edit the laser range observations. 

For each observing station the quality of the interpolated ranges as 
reflected through the qualitative pattern of the Chebychev residuals is 
relatively homogeneous for all of the arcs recorded during the MERIT 
Main Campaign. Thus, in order to get a feeling for the quality of the 
interpolated ranges we have plotted in Figs. 14 and 15 the Chebychev 
residuals for those arcs of Figs. 8 and 13 which have been marked with 
an asterisk next to their numeric station ID’s. If a specific station has 
observed more than one of the arcs shown in Figs. 8 and 13, then the 
arc whose Chebychev residuals are shown in Fig. 14 or 15 is also 
marked just above its starting point with an asterisk. Table 5 lists, for 
the same arcs shown in Figs. 14 and 15, the RMS values of the 
Chebychev residuals together with the number of the available 
observations before and after the data-snooping procedure was applied. 
The condition numbers also shown in this table refer to the normal 
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Fig. 13 Arc overlap of three Lageos passes (American 
stations). 
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Fig, 13 (conVd) 

equations matrix that was used to estimate the coefficients of Iho 
least-squares Chebychev polynomial. The condition numbers missing 
from this table correspond to the stations whose actual ranges were 
used to generate both the SR and SRD observables. It is evident from 
this table that for some stations rejection of less than 10% of their 

observations reduces the RMS of the Chebychev residuals from 
unacceptable levels down to the expected accuracy of the observations. 
Furthermore, the large condition numbers associated with stations 7210 
and 7121 are caused by the strong irregularities in the distribution of 
the data gaps. These strong irregiilarities usually occur when 

malfunction of the laser instrument has resulted in concentrated 
erroneous observations which in turn are rejected by the data-snooping 
procedure. For instance, the observations of stations 7210 and 7121 for 
the arcs shown in Fig. 15 are affected by four large gaps and one large 
gap respectively. The word large is used here in the sense that 

although the data gaps are usually shorter than 60 seconds they are 
still large as compared to the other gaps affecting the remaining 

observations. 
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Table 4, Monthly F'recision Estimates (cm) for the American 
Stations (Chebychev Fi tting/ After Data Snooping) 
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Chebychev residuals for stations 7907| 7105, 7086, 7109, 7110, 7886 (after data 
snooping)* 





ITMD 8M082M ITMO 840824 

OBSER. 296 RMS =0.06 OBSER. 9474 RMS 
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Fig. 14 (cont'd) 




Fig. 14 (cont'd) 




ITMD 840407." ITMD 840125 

OBSER. 107 RMS =0.11 OBSER. 811 



Pig. 15 Chebychev residuals for stations 7112, 7265, 7121, 7210, 7220, 7122, 7062 (after data 
snooping). 
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Fig. 15 (cont’d) 

The large condition numbers, caused by the strong irregularities in 
the distribution of the gaps, might deteriorate the precision of the 
interpolated ranges since the number of the significant digits lost in the 
inversion of the normals is approximately equal to the base 10 logarithm 
of the condition number (Forsythe et aL, 1967). Even if the numerical 

instability is not bad enough to affect the quality of the interpolated 
ranges there might still exist large approximation errors in the 
neighborhood of the large gaps. Therefore, it was decided to either 
reject those arcs or to break them down into several subarcs with even 
distribution of gaps (i.e., small condition number) provided that for each 
subarc there are enough available observations to be effectively 
interpolated. 

Sufficient representation of a relatively large number of 
observations requires a high-degree Chebychev polynomial. 
Determination of its coefficients through a least-squares adjustment does 
not cause any numerical instability, as might have been expected, 
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Table 5. Precision of Chebychev Interpolation 
(Before and After Data Snooping) . 


Station 

IDs 

RMS 

before 

D.S.*(m) 

No. of obs. 
before 

D.S. 

RMS 
after 
D.S. (m) 

No. of obs. 
after 
D.S. 

Condition 

Number 

Reference 

fibres 

7907 

26.9 

648 

0.13 

620 

600 

8, 

14 

7105 

7.35 

5314 

0.02 

5116 

561 

8, 

14 

7086 

0.07 

297 

0.06 

296 

— 

8, 

14 

7109 

0.04 

9611 

0,02 

9474 

31 

8, 

14 

7886 

0.07 

2720 

0.07 

2720 

102 

8, 

14 

7110 

0.23 

6681 

0.02 

6570 

34 

8, 

14 

7122 

0.16 

2088 

0.15 

2072 

48 

13, 

15 

7220 

0.09 

686 

0.09 

685 

106 

13, 

15 

7062 

0.70 

613 

0.09 

543 

— 

13, 

15 

7265 

0.09 

812 

0.09 

811 

— 

13, 

15 

7112 

23.18 

136 

0.11 

107 

— 

13, 

15 

7210 

0.07 

185 

0.05 

167 

4249 

13, 

15 

7121 

15.88 

1694 

i 

0.14 

1653 

1058 

13, 

15 


♦Data snooping 


because of the near orthogonality properties of the Chebychev 
polynomials (see Section 3.5.2). For instance, to sufficiently represent 
9474 observations of station 7109, a 22— degree Chebychev polynomial is 
required. In this case the condition number of the normals has a very 
small value (i.e., 31, from Table 5), and therefore numerical instability 
and ill-conditioning associated with large condition numbers (ibid.) are 
not a concern in spite of the fact that a high— degree Chebychev 
polynomial was used. 

All the information acquired in regard to the interpolation of the 
laser ranges was implemented in the locally developed software that was 
subsequently used to generate the SR and SRD observables. With this 
software the passes containing arcs coobserved by two or more stations 
were identified and isolated together with the starting and ending 
epochs of each of the arcs involved. More specifically, for the MERIT 
Main Campaign there have been identified 536 such passes observed by 
the American stations shown in Fig. 12. For each of these passes the 
observations recorded within an arc by all of the stations involved were 
edited by using the data-snooping procedure and by taking into 
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consideration all of the necessary precautions in regard to the 
magnitude, the distribution of the gaps, and the resulting numerical 
instability* However, effective editing of the laser ranges requires 
knowledge of the degree of the Chebychev polynomial that would 
sufficiently represent the laser ranges recorded within any of the arcs 
involved* This is accomplished with the help of a statistical test which 
is explicitly outlined in Section 3*6.2. According to this test, any 
adjustment for which the hypothesis Ho is rejected leads to another 
adjustment (see Section 3.6.2). To eliminate these additional 
adjustments, and therefore to reduce the computing time, it was decided 
to determine a priori the degree of Chebychev polynomials that 
sufficiently represent the observations of any of the stations involved 
in our study (see Fig. 12). Accordingly, if the number of available 
observations recorded by a specific station falls into a certain window 
then for that station an appropriate degree of the Chebychev polynomial 
is assigned (see Table 6). These degrees are predetermined for 
specified windows and for all of the stations involved on the basis of 
the actually observed ranges. Table 6 lists for prespecified windows 
and for all of the stations involved the degree of the Chebychev 
polynomials that sufficiently represent the ranges falling into any of 
those windows. 

Next the software proceeds with the identification of the passes 
containing arcs that have been coobserved by four or more stations. 
For each of those passes the exact overlapping times of the arcs 
involved are identified and subsequently Simultaneous Ranges (SR) are 
generated along the guidelines described in Section 3.6. For instance, 
stations 7110, 7122 7220 and 7062 have arcs with quasi-simultaneous 
observations (see Fig. 16). It is obvious from Fig. 16 that 

quasi-simultaneity for all four stations occurs only in the intervals (be) 
and (hi). The epochs b and c designate the starting and ending epochs 
of the first arc of station 7062, while the epochs h and i designate the 
starting epoch of the second arc and the ending epoch of the fourth 
arc of stations 7220 and 7062 respectively. Since station 7062 recorded 
the least number of observations for both arcs be and hi respectively, 
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Table 6. Degree of Chebychev Polynomials. 


N. No. of 

obs. 

StatioiNv 
IDs X 

0 - 
500 

500 - 
1000 

1000 - 
2000 


3000 - 
5000 

5000 - 
9000 

> 9000 

7907 

16 

17 

17 

18 

19 

21 

22 

7105 

12 

16 

16 

18 

19 

21 

22 

7086 

12 

17 

17 

18 

19 

21 

22 

7109 

12 

17 

17 

18 

19 

21 

22 

7886 

12 

13 

15 

17 

17 

21 

22 

7110 

12 

16 

16 

19 

20 

21 

22 

7122 

14 

19 

19 

20 

21 

21 

22 

7220 

11 

11 

14 

14 

16 

21 

22 

7062 

16 

18 

18 

19 

19 

21 

22 

7265 

11 

13 

15 

18 

19 

21 

22 

7112 

14 

19 

19 

20 

21 

21 

22 

7210 

11 

19 

19 

20 

21 

21 

22 

7121 

14 

19 

19 

20 

21 

21 

22 


the four station events (i.e., exact simultaneous observations from four 
stations) for these two arcs were generated by interpolating all the 
available observations from the corresponding arcs of stations 7110, 7122 
and 7220 respectively. Although in Fig. 16 we have only four-station 
arc overlaps, there may exist arc overlaps including five, six or even 
seven stations. In such cases all possible seven-, six-, five- or 
four-station events are generated, provided that any duplication that 
may occur is to be avoided. 
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Fig. 16 Arc overlap of a Lageos pass involving four 
American stations. 

Finally all the passes with quasi-simultaneous observations from two 
or more stations were employed to generate the 3RD observables. For 
instance, the station pair 7122-7110 shown in Fig. 16 has 

quasi-simultaneous observations for the intervals (ad), (ef) and (gj) 
respectively (see Fig. 16). For the first interval (ad) the 3RD 
observables are generated by interpolating the observations of the first 
arc of station 7110 at the observing epochs of station 7122, while for 
the intervals (ef) and (gj) the 3RD observables are generated by 
interpolating all the available observations of station 7122. Thus, having 
all possible four-, five-, six- and seven-station events together with all 
possible SRD's, the baselines are estimated through the geometric and 
3RD methods respectively. This is the subject of the next chapter. 

118 

viIGlNAL PAGE JS 
OY POOR QUALITH 



Chapter 4 


BASELINE ESTIMATION 

4.1 BASELINE ESTIMABILITY 

Baseline estimability in connection with satellite geodesy is very closely 
related to the concept of estimable parameters. In fact, if optimum 
geometric configurations are fulfilled and enough observations are 
available, the baselines can be estimated with a precision compatible to 
that of the observations. This is possible only because the baselines 
form a set of estimable parameters. This is true for both the range 
geometric mode and the 3RD semidynamic/dynamic mode methods. In 
general, baselines estimated through either geometric or dynamic mode 
methods are estimable only if through the appropriate adoption of the 
necessary constants and units the scale is implied by the measurement 
system being employed. In such cases the baselines estimated with a 
dynamic mode method are referred to as best estimable parameters. The 
use of this term is justified since out of all possible estimable 
parameters associated with either semidynamic or dynamic mode methods, 
the baselines are recovered with substantially reduced a posteriori 
variances (i.e., an order of magnitude) (Van Gelder, 1978). 

4.2 STEADY STATE RESPONSE OP THE GEOMETRIC AND SRD METHODS 
The nature of as well as the spatial and temporal distribution of the 
available observations quantify the inherently present information for 
each of the estimable parameters involved (Sections 2.2, 4.3 and 4.4). 

The process of recovering estimable quantities from any given set of 
observations is referred to as an inversion process. This process is 
effectively realized via an estimation method such as a least squares 
adjustment (Sections 2.1.2 and 2.2.6). The available set of observations 
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forms the input to the inversion process, while the recovered estimable 
quantities constitute the response (i.e*, the output) of this process to 
the given set of observations. 

Different inversion processes have different responses to a specific 
set of observations and adjusted quantities. Nevertheless, as far as 
satellite geodesy is concerned, these responses cannot be meaningfully 
differentiated in the light of the current observational distribution and 
accuracy. 

The steady state response of an estimation method has been reached 
if extension of its input with additional observations does not contribute 
any additional information to the estimable quantities being recovered. 
In the present study, steady state responses of both the geometric and 
the SRD methods is sought through the extension of their input with 
more and more observations. This process continues up to that point 
where additional observations do not affect the recovered baselines 
beyond the level of accuracy implied by the sophistication of the models 
employed and by the accuracy of the available observations. With iho 
assumption that the accuracy of the baselines, recovered from laser 
range observations collected during the Main MERIT Campaign, cannot 
exceed the 1 cm level, the steady state response of both the SRD and 
the geometric methods is said to have been reached if additional 
observations do not affect the recovered baselines at the 1 cm accuracy 
level. Furthermore, we refer to steady state response of the SRD 
method” if the input observables are the dynamically modeled SRD 
observables, while we refer to ” steady state response of the geometric 
method” if the input observables are the geometrically modeled 
simultaneous ranges (see Chapter 2). 

In many circumstances the steady state response of either the 
geometric or the SRD methods may not be possible because their 
response either diverges or oscillates. Such a response, however, can 
be reached if the geometric and/or physical characteristics of either the 
geometric or the SRD methods are changed accordingly. 
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In the geometric method, for instance, near singularity cases 
(Section 2.1.3) result in an extremely ill-conditioned normal equation 
matrix which in turn leads to a divergent response. Divergence is 
reversed to convergence if the input of the geometric method is 
extended to include more and more observations collected by stations 
located well away from the critical surfaces and the critical curves 
(Blaha, 1971). However, if the distribution and number of available 
observations is not good enough to warrant a steady state response, 
then such a response can be reached via the implementation of 
appropriately chosen constraints. This is accomplished either by 
constraining baselines in the geometric method (Section 4.3.2) or by 
increasing the lengths of continuous integration in the SRD method 
(Section 4.4.2). 

Increasing the length of continuous integration in the implementation 
of the SRD method results in a "faster” steady state response in regard 
to baseline estimation. The term "faster" is used to indicate that the 
steady state response of the SRD method is reached on the basis of a 
substantially reduced number of observations compared to those needed 
to reach this steady state response via the short arc solutions (Section 
4.4). This is the result of the geometric strength implied by the long 
integration periods and manifested in the reduced order of the normal 
equation matrix and in the larger values of its corresponding diagonal 
elements. 

Long integration periods constitute a potential source of 
"instability" of the normal equation matrix. The term "instability" is 
used to express the existence of high correlations among the adjusted 
parameters. High correlations (greater than .99) raise the warning flag 
of ill-conditioning of the normal equation matrix and a possible 
divergent response of the SRD method. Employing long integration 
periods independently of what the baseline pass geometry is, the state 
vectors of the passes that are days away from the corresponding 
adjusted initial state vector are largely insensitive to that initial stale 
vector. Therefore, for those passes there exists a tendency for litKiar 
dependence among the corresponding columms of the state transition 
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matrix Y*(t) (equation 2-63), This tendency further contributes to the 
instability of the normal equation matrix. The instability of the normal 
equation matrix gets even worse when single long baselines (longer than 
about 4000 km) are estimated via the SRD method. For those baselines 
there is a tendency for th€) coobserved passes to concentrate in the 
areas between the end points of the baseline (Section 4.4.1). 

With the weather-dependent laser range observations it is very 
likely that the stations consitituting the end points of the estimated 
baseline have coobserved only one pass falling within a specific 
integration period. Under these circumstances the length of continuous 
integration will be reduced to that of the duration of one single pass 
which lasts about one hour for the Lageos satellite. Such short 
integration lengths implemented in the same solution with integration 
lengths of up to one week (Section 4.4) constitute a potential source of 
instability because of the resulting inhomogeneity in the structure of 
the normal equation matrix. In the design matrix, this inhomogeneity 
manifests itself by the many zero entries affecting the columns 
associated with the initial state vectors of those short arcs. Therefore, 
care should be exercised to avoid extreme circumstances that may result 
in an algorithmically singular normal equation matrix (Section 4,4.2). 
This can be avoided either by employing homogeneous integration 
lengths and/or by incorporating in the same SRD solution observations 
from many stations. This would improve the) stability of the normals 
because of the strength implied by the geometry of the additional 
observations. This, however, is true only if four or more stations are 
involved since as it was mentioned in Section 2.1,3, at least four stations 
and six targets are generally necessary for a nonsingular range 
geometric solution. Therefore, in the single baseline solutions the 
stability of the normals is mainly controlled by the constraint imposed 
on the satellite to move along a six-parameter orbit (Chapter 2). 

In the present study only single baseline solutions via the SRD 
method have been performed. This was decided upon because, for 
moderate baseline lengths (<3500 km), the single SRD baseline sohations 
seem to be simple, fast and very accurate. 
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The maximum length of continuous integration in the implementation 
of the SRD method should be chosen not only to assure the steady state 
response of the SRD method but also to warranty minimum propagation 
of the accumulated residual orbital biases into the recovered baselines. 
These biases are the ones not eliminated by the nature of the SRD 
observables which reduces substantially the accumulated orbital biases 
with an almost total cancellation of the accumulated radial biases (Pavlis, 
1982). Consequently, the choice of the integration lengths is limited 
only to those which result in an a posteriori standard deviation of unit 
weight the value of which is close to unity. Thus, in seeking the 
steady state response of the SRD method we start out with short arc 
solutions. If with short arc solutions the available observations are not 
enough to warranty a steady state response, the integration lengths are 
steadily increased up to those resulting in a steady state response. If 
the a posteriori standard deviation of unit weight approaches unity 
before such a response has been reached, then the corresponding 
baseline is not estimated. However, any integration lengths leading to a 
steady state response will be adopted even if the a posteriori standard 
deviation of unit weight is smaller than one. This simply means that the 
baselines via the SRD method are estimated on the basis of the minimally 
required integration lengths. 

The effectiveness of the a posteriori standard deviation of unit 
weight to control the maximum length of continuous integation in the 
single SRD baseline estimation depends heavily on the assumption of 
having weighted the observations properly. This assumption has been 
elaborated on and it has also been well justified in Chapter 3. 

4.3 BASELINE ESTIMATION VIA THE STEADY STATE RESPONSE OP THE 
GEOMETRIC METHOD 

In the absence of ill-conditioning that may result from nearly critical 
configurations, the steady state response of the geometric method will 
be reached if a large enough number of observations is available. This 
is a consequence of the fact that the weight coefficient matrix (Qx) of 
the adjusted parameter vector "improves" as additional observations are 
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incorporated into the geometric solution (Blaha, 1971). The term 
"improves" is used to indicate that with the same adjusted parameter 
vector the matrix (Q$ - Qx^) is a positive (semi) definite matrix where 
(Qx) and (Qx^) denote the weight coefficient matrices of the adjusted 
parameter vector before and after the inclusion of additional 
observations. Since (Qx - Qx^) is a positive (semi) definite matrix, it 
follows that (ibid.) 

Tr (Qj^) < Tr (Qg) (4~1) 

indicating that the precision of the adjusted parameters improves only if 
the additional observations are drawn from the same population, ther‘eby 
making it possible to assume that the a priori variance of unit weight is 
the same for both solutions and it is equal to one. Since the 

mathematical model employed in the minimum constraint geometric 
adjustments is an almost error free model, the a posteriori variance of 
unit weight tends towards one at the steady state. If these a posteriori 
variances of unit weight are very close to unity so they preserve the 
relative structure of the weight coefficient matrices of eq. (1-1), thf> a 
posteriori variance of the adjusted parameters decreases. In the 
geometric mode adjustments, the adjusted parameters are the Cartesian 
coordinates of the ground stations. 

In the present study we are particularly interested in estimating 
baseline lengths and therefore the a posteriori variances of the 
Cartesian coordinates have been mapped bac:k into the estimated 
baselines. As the a posteriori variance of the estimated Cartesian 
coordinates improves with the incorporation of additional observations, 
the a posteriori variance of the estimated baselines will also improve 
only if the nonlinearity of the model employed allow suc:h an 
improvement to take place. Therefore, by including additional 
observations into the geometric solutions a point will be reached where 
the estimated baseline length will not change beyond the 1 cm level. 

It turns out that before reaching the steady state, the a posteriori 
standard deviations of the estimated baselines do not reflect the change 
of their length as additional observations are incorporated into the 
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solution. For instance, it is possible for an estimated baseline, the a 
posteriori standard deviation of which is 5 cm, to change its length by 
as much as 76 cm with the incorporation of additional observations 
(Table 8, baseline 7110-7122). This means that before reaching the 
steady state response, the a posteriori standard deviations of the 
estimated baselines do not reflect their accuracy but rather they 
indicate how far away the solution is from its steady state. If, however, 
the steady state response has been reached, the corresponding a 
posteriori standard deviations assume millimeter-level values (Table 8). 
Therefore, in the present study the a posteriori standard deviation of 
the estimated baselines are used only as indicators showing whether or 
not the steady state response has been reached. Successful utilization 
of these indicators assumes close to unity a posteriori variances of unit 
weight so they preserve the relative structure of the weight coefficient 
matrices of the recovered parameters as additional observations are 
incorporated into the solution (eq. 4-1). 

If steady state has been reached, it is furthermore assumed that the 
corresponding baseline has been determined at the 1 cm accuracy level. 
This is well Justified only if the following three assumptions are 
fulfilled: (1) the motion of the observing stations has been properly 

modeled for the time span of the observations, (2) the steady state 
response implies strong geometry, and (3) the accuracy of the available 
observations allows recovery of the baselines at the 1 cm level (see 
Section 4.2). As of the geometric modeling itself, the only errors 
affecting it arise from using the three-dimensional Euclidean space 
formulation rather than the four-dimensional post-Newtonian formalism 
(Moyer, 1971; Bjerhammar, 1985) which of course does not affect the 
estimated baselines at the 1 cm level (Moritz, 1979). 

4.3.1 Geometric Strength of the Available Observations 
In the geometric approach the observed satellite positions are treated as 
auxiliary independent points in space (see Section 2.1); therefore the 
strength of any minimum constraint geometric solution depends entirely 
on the geometric strength implied by both the amount and the 
distribution of the available observations. Thus, any meaningful 
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interpretation of the results obtained in geometric mode adjustments 
should be preceded by an analysis aiming to assess the strength of the 
geometry involved in each of those adjustments. 

The analysis of the geometry involved in the geometric adjustments, 
presented later in this section, is based on the examination of Table 7. 
The first column of this table lists the ID of the stations involved, the 
second, third and fourth columns contain the number of observations 
per station when all the simultaneous events from Sept. 1983 to May 
1984, from Sept. 1983 to Aug. 1984, and from Sept. 1983 to Oct, 1984 are 
considered respectively. The bottom part of the second and third 
columns contain the total number of observations, the degrees of 
freedom and the a posteriori variance of unit weight obtained from 
minimum constraint solutions on the basis of the data shown in the 
corresponding column. The bottom part of the third column contains the 
same information for two types of solutions, one minimum constraint and 
one overconstraint obtained using the data listed in that column. 

The geometric strength in each of the above solutions is primarily 
drawn from the stations having the most observations. In assessing 
this strength we assume that stations 7886 and 7220 coincide with 

stations 7109 and 7110 because they are 8 m and 15 m away from each 
other. With this in mind it can be easily deduced from Table 7 that 86% 
of the available observations has been recorded by stations 7105, 7109, 
7110 and 7122. A geometric solution, the strength of which is primarily 

drawn from four stations, tends to be sensitive to how close these 

stations are from their best fitting plane (Blaha, 1971) because with 

stations close to forming a plane six are needed for a nonsingular space 
range network (see Section 2.1.3). Thus, when six or more stations are 
involved the solutions are not sensitive to the closeness of those 
stations from the best fitting plane. 

In the geometric solutions shown in Table 7 the additional 11% of 
the available observations has been recorded by stations 7112, 7086 and 
7265. Therefore, 97% of the observations has been recorded by seven 
stations. The existence of these three additional stations reduces the 
sensitivity of the solutions to how close stations 7105, 7109, 7110 and 
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Table 7 Global Statistics of the Geometric Adjustments 


Station 

ID 

Sep 83 - May 84 

No. of Observations^^) 

Sep 83 - Oct 84 

Sep 83 - Aug 84 

7105 

7143 

18990 

19214 

7109 

10198 

22784 

23936 

7112 

2884 

3467 

3467 

7122 

10996 

11284 

12212 

7220 

1969 

1969 

1969 

7110 

11549 

24200 

25352 

7062 

841 

841 

841 

7086 

1412 

4400 

4624 

7907 

176 

245 

245 

7082 

299 

299 

299 

7210 

712 

756 

1684 

7265 

4395 

4395 

4395 

7886 

— 

11859 

11859 


N = 52574 

N = 105489 

N = 110097 


DF = 14519 

DF = 29478 

DF = 30630 


= 1.18 

<To* = 1.03 

<ro* = 1.03 




DF = 3063l('“) 




= 1.05(=^) 


Minimum constraint solution, coordinates fixed(^); X,Y,Z for 
7109; X,Y for 7122; Z for 7105. 

Overconstraint solution, coordinates fixed(^): X,Y,Z for 7109; 
X,Y,Z for 7122; Z for 7105. 

Coordinates fixed to those of (CSR)85L01 (Section 4.4.2) 

N = total nvunber of observations 
DF = degrees of freedom 

<To^ = a posteriori variance of unit weight 


7122 are from their best fitting plane. The remaining four stations, 
7062, 7907, 7082 and 7210, have contributed only the last 3% of the 
available observations. Therefore, their contribution to the geometric 
strength is very minor as compared to that of the previous stations. 

Next we assume that stations 7112, 7086 and 7265, which contribute 
the 11% of the available observations, make the geometric adjustments 
insensitive to how close stations 7105, 7109, 7110 and 7122 are from their 
best fitting plane. Even if this were true, the strength of the 
geometric solutions tends to still be weak because the stations that have 
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Table 8 Baseline Steady State Response of the Geometric Solution 


Baseline 

No. of 
Observ. 

Correlation 
> 0.80 

Length (m) 

Data 
Set ‘ 

Solution 
Type 2 

7109-7110 

9,186 

None 

883601.637 ±0.02 

A 

1 


21,772 

ft 

.608 * 0.02 

B 

1 


22,924 

ft 

.661 ± 0.02 

C 

1 


tf 

ft 

883601.661 ‘0.02 

C 

2 


»» 

ff 

883602.245 ‘ 0.009 

C 

3 

7109-7265 

3,363 

ft 

627043.412 ‘0.02 

A 

1 


ff 

M 

.452 ‘ 0.01 

B 

1 


M 

ff 

.535 ‘ 0.01 

C 

1 


tf 

ff 

.535 ‘ 0.01 

C 

2 


tf 

ff 

.988 ‘ 0.005 

c 

3 

7109-7886 

11,859 

None 

7.746 ‘ 0.002 

B 

1 


ft 

ff 

.746 

C 

1 


ft 

ff 

. 746 

C 

2 


ft 

ft 

.746 

c 

3 

7109-7122 

8,644 

None 

2280712.335 ‘0.07 

A 

] 


8,932 

ft 

2.700 ‘ 0.05 

B 

1 


9,860 

ft 

3.188 ‘ 0.05 

C 

1 


ft 

ft 

3.188 ‘ 0.05 

C 

2 


ff 

ft 

4.949 ‘ 0.0005 

C 

3 

7110-7122 

10,060 

PZfZs=0.998 

1437137.428 ‘0.05 

A 

1 


10,348 

ft 

.780 ‘ 0.04 

B 

1 


11,276 

If 

8.187 ‘ 0.03 

C 

1 


It 

^XpXg— 0, 996 

8.187 ‘ 0.03 

C 

2 


tf 

None 

9.288 ‘ 0.009 

C 

3 

7110-7220 

1,576 

/^XpXs=0, 801 

15.225 ‘ 0.006 

A 

1 


PYfYs=0. 990 
PZfZs=0.998 


ff 

a>YfYs=0.987 

/^ZfZs=0.998 

.221 

± 

0.005 

B 

1 

ff 

PYfYs=0.981 

PZfZs=0.997 

.218 

± 

0.005 

C 

1 

ff 

PXpXs=0.997 
PyfYs=0. 975 

PZfZs=0. 800 

.218 

± 

0.005 

C 

2 

ff 

PYfYs=0.965 

.208 

± 

0.005 

C 

3 

7110-7265 3,866 

PYfYs=0. 91 
PZfZs=0.996 

274069.453 

± 

0.01 

A 

1 

ff 

PyfYs-0.887 

PZfZs=0.995 

.383 

± 

0.008 

B 

1 

ff 

PyfYs=0. 850 
PzfZs=0.994 

.355 

± 

0.008 

C 

1 

tf 

PxfXs=0.994 

.355 

± 

0.008 

C 

2 

ff 

None 

.474 

± 

0.007 

C 

3 
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Table 8 (cont’d) 


Baseline 

No. of 
Observ . 

Correlation 
> 0.80 

Length 

(m) 

Data 
Set ‘ 

Solution 
Type * 

7110-7886 

11,859 

None 

883605.698 

± 

0.02 

B 

1 


tf 

tf 

.751 

± 

0.02 

C 

1 


ft 

tf 

.751 

± 

0.02 

C 

2 


ft 

ff 

6.335 

± 

0.009 

C 

3 

7122-7265 

4,184 

PzfZs=0.995 

1663980.848 

± 

0.05 

A 

1 


fl 

=0.994 

1.161 

± 

0.04 

B 

1 


tf 

=0.993 

1.555 

t 

0.04 

C 

1 


ft 

^XpXs=0.993 

1.555 


ft 

C 

2 


ft 

None 

2.823 

± 

0.005 

C 

3 

7122-7886 

0 

None 

2280718.021 

± 

0.05 

B 

1 


0 

ft 

.509 

t 

0.05 

C 

1 


0 

tf 

18.509 

± 

0.05 

C 

2 


0 

ft 

20.269 

± 

0.002 

C 

3 

7220-7265 

0 

PyfYs=0.90 

274066.158 

t 

0.010 

A 

1 



PzfZs=0.994 







0 

PyfYs=0.874 

.090 

± 

0.009 

B 

1 



/^ZfZs=0.993 







0 

PyfYs=0.833 

.064 

± 

0.008 

C 

1 



PzfZs=0.991 







0 

PxfXs=0.991 

.064 

± 

0.008 

C 

2 


0 

None 

.189 

t 

0.007 

C 

3 

7265-7886 

0 

None 

627048.351 

± 

0.01 

B 

1 


0 

ft 

.434 

t 

0.01 

C 

1 


0 

tf 

.434 

± 

0.01 

C 

2 


0 

ft 

.887 

± 

0.006 

C 

3 


’ Data Sets: A Sep 83 - May 84 
B Sep 83 - Aug 84 
C Sep 83 - Oct 84 


* Solution Type: 

1 Minimum constraint solution, Cartesian coordinates fixed: 
X,Y,Z for 7109; X,Y for 7122; Z for 7105 

2 Minimum constraint solution, Cartesian coordinates fixed: 
X,Y,Z for 7109; Y,Z for 7122; Z for 7105 

3 Overconstraint solution, Cartesian coordinates fixed: 
X,Y,Z for 7109 ^ 7122; Z for 7105 


recorded the 91% of the observations (i.e., 7105, 7109, 7122, 7220, 7110, 
7062, 7265 and 7886) are concentrated around two intersecting lines 
defined by station 7109 with 7122 and 7122 with 7105 (see Pig. 8). Since 
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two intersecting lines belong in the family of second-order curves, it is 
reasonable to expect that near singularity B (Section 2.1.3) tends to 
weaken the strength of those solutions. If we furthermore assume that 
the 9% of the observations recorded by stations 7082, 7086, 7112, 7907 
and 7210 which are located well away from these two intersecting lines, 
make the geometric solutions insensitive to the previously mentioned 
near singularity B, it is still possible that near singlarity C might be 
present. This is a consequence of employing a network of relatively 
large extent. With such a network the simultaneously observed satellite 
positions tend to concentrate on the area extended above the middle 
part of the network and therefore to be closer to a plane. This in turn 
would lead to near singularity C because off-plane targets are needed to 
avoid this type of singularity (Section 2.1.3). 

It is evident from the above discussion that the geometric strength 
implied by the available observations is relatively weak. Therefore the 
minimum constraint geometric solutions will be strongly influenced by 
the distribution of the available observations thereby making it difficult 
to recover the relative geometry of the ground stations (next section). 

4.3.2 Baseline Results 

On the basis of the data listed in the last three columns of Table 7, we 
have performed five least squares geometric mode adjustments. Four of 
these adjustments are based on minimum constraints, and the fifth is 
based on constraining one additional Cartesian coordinate than those 
required for a minimum constraint solution. More specifically, the first 
three minimum constraint geometric adjustments have been obtained by 
using all the events from Sept. 1983 to May 1984, from Sept. 1983 to 
Aug. 1984 and from Sept. 1983 to Oct. 1984. An event occurs when four 
or more stations observe the same satellite position. In these three 
adjustments, the same mimimum constraints have been used as they are 
implied by fixing the X,Y,Z Cartesian coordinates of station 7109, the X 
and Y coordinates of station 7122, and the Z coordinate of station 7105. 
The fourth adjustment was obtained by using all the available events 
and by fixing the Y and Z rather than X and Y coordinates of station 
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7122. This adjustment was necessary to confirm further the weakness 
of the geometry that seriously affects the minimum constraint geometric 
solutions. Finally, the overconstraint solution performed on the basis of 
the data listed in the fourth column of Table 7 was necessary to reach 
the steady state response for longer baselines. The results for all of 
the adjustments decribed above are listed in Tables 7 and 8. Table 8 
contains for each baseline and for all of the solutions performed the 
number of observations per baseline, the correlations that are greater 
than or equal to 0.80 and the estimated baseline lengths followed by 
their a posteriori standard deviations. Out of all of the estimated 
baselines, only those are listed for which the a posteriori standard 
deviations in the overconstraint solution have reached the millimeter 
level. The symbol PZpZs* in Table 8, designates the correlation between 
the Z coordinate of the first (F) station and the Z coordinate of the 
second (S) station of the corresponding baseline. 

A close inspection of the minimum constraint solutions clearly 
reveals that the estimated baselines whose a posteriori standard 
deviation is at the centimeter level change their length by as much as 
70 cm to 80 cm when additional observations are incorporated into the 
solutions (baselines 7110-7122 and 7122-7265). The number of additional 
observations incorporated in each solution as compared to the previous 
one is easily deduced from the information given in Table 7 for each one 
of those solutions. 

The baselines whose length is smaller than 1000 km change their 
length to within -10 cm to 12 cm with a tendency of positive increase. 
This positive increase, which is clearly pronounced for the longer 
baselines (7110-7122 and 7122-7265) takes place toward the correct 
length of the corresponding baselines. This is easily seen by comparing 
these lengths with the ones obtained in the overconstraint adjustment 
(i.e., Solution C 3). The fact that the lengths of the baselines, which 
have been estimated at millimeter precision level via the overconstraint 
adjustment are very close to their true lengths is elaborated on at the 
end of this section, and it has also been confirmed by comparing their 
lengths with those obtained via the SRD method as well as with those 
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computed by other computational centers such as CSR and ZIPE via 
dynamic long-arc solutions (see next section). 

The geometry implied by the number and the distribution of the 
available observations manifests itself in the minimum constraint 
solutions shown in Table 8. This geometry is not strong enough to 
warrant a steady state response since all the baselines but 7110-7220 
and 7109-7886 change their length by several centimeters with the 
incorporation of additional observations. The existence of weak 
geometry is the result of the expected near singularities described in 
the previous section. This weak geometry is also confirmed by the high 
correlations prevailing among the station coordinates that were not 
constrained in the implementation of the minimum constraint solutions. 
This weakness is further confirmed by using all the available events 
from Sept. 1983 to Oct. 1984 and by changing the minimum constraints 
from (1) to (2) (Table 8). In solution (2) the correlations prevail among 
the X coordinates of the stations for which the correlations in solution 

(1) were high among their Z coordinates. It is interesting to note that 
the correlations among the X coordinates in solution (2) are almost the 
same as those among the corresponding Z coordinates of solution (1). 
This is the result of fixing the X and Y coordinates of station 7122 in 
solution (1) and the Y and Z coordinates of the same station in solution 

(2) . The reduction of the correlations among the Y coordinates of the 
corresponding stations from (1) to (2) simply reveals that the minimum 
constraints (2) have better stability characteristics than those of 
minimum constraint (1). 

The correlations among the coordinates of station 7110 and 7220 are 
also high because these stations are very close to each other, about 15 
m apart. This, in light of weak geometry, makes their separation 
difficult thereby bringing the corresponding correlations to high levels. 
The correlations among the coordinates of station 7886 and the 
coordinates of all of the other stations involved are small because 7886 
is only about 8 m away from 7109 which was held fixed in the 
implementation of both solutions (1) and (2). The estimated baseline 
lengths and their a posteriori standard deviations remained of course 
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the same for both solutions because the baselines are invariant 
quantities under any minimum consraints and as such their estimated 
lengths together with their variances should not change from one 
minimum constraint solution to the next. 

The above discussion evidently reveals that the geometry implied by 
all of the simultaneously observed satellite positions is weak thereby 
making the minimum constraint solutions susceptible to the distribution 
of the available observations. For reasons mentioned in the previous 
section, the simultaneously observed satellite positions tend to 
concentrate on the area extended above the middle part of the employed 
network (see Fig. 12). This observational coverage together with the 
weak geometry will imply, via the minimum constraint solutions, a range 
space network with a tendency to shrink towards its center and more 
specifically towards the area where most of the available observations 
are concentrated. This simply means that the scale of the recovered 
range space network is not properly implied by the geometric strength 
of the available observations. This fact is also confirmed by the 
positive increase of the longer baselines when additional observations 
are incorporated into the minimum constraint solutions. Implementation 
of the scale in the geometric solution has been attempted by fixing the 
third coordinate of station 7122 (solution (3)) in addition to the 
coordinates fixed in the implementation of the minimum constraint 
solutions. By fixing the third coordinate of this station we implicitly fix 
the length of baseline 7109-7122. 

As mentioned earlier in this section. Table 8 contains the results 
only for those baselines for which the a posteriori standard deviations 
in the overconstrained solution are smaller than 1 cm. Although all of 
the remaining baselines have been estimated in all solutions shown in 
this table, they are not included there because the steady state 
response for those baselines is assumed not to have been reached. The 
reason for claiming this is based on the fact that baselines 7110-7220 
and 7109-7886 are the only ones that do not change their lengths at the 
centimeter level either when additional observations are incorporated 
into the solution or when additional constraints are applied. 
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Furthermore, only these two baselines have been estimated in the 
minimum constraint solutions with an a posteriori standard deviation at 
the milimeter level thereby suggesting that a steady state response in 
the geometric mode environment assumes millimeter level a posteriori 
standard deviations of the estimated baselines. 

Application of an additional constraint brings the a posteriori 
standard deviation of all the remaining baselines shown in Table 8 to 
the millimeter level thereby implying that for those baselines the steady 
state response has been reached via the implementation of the scale in 
the resulting range space network. 

This was expected because in the geometric solutions the baselines 
constitute a set of estimable quantities and as such constraining one of 
them (i.e., 7109-7122) will lead, apart from nonlinear terms, to reduced 
unsealed standard deviations for the remaining baselines. Their scaled 
standard deviations (i.e., a posteriori standard deviations) will also be 
reduced if the scaling factors (i.e., the a posteriori standard deviations 
of unit weight) are close enough to unity so they preserve the relative 
structure of the corresponding unsealed standard deviations. In fact 
the a posteriori standard deviation of unit weight in the overconstrained 
solutions is 1.05, just slightly larger than that of the corresponding 
minimum constraint which is 1.03, making it possible to preserve in 
those solutions the same relative structure of the scaled variance- 
covariance matrices as that of their corresponding weight coefficient 
matrices (Table 8). This simply means that the additional constraint did 
not distort the geometry but rather made it stronger, and therefore it 
was possible to reach through the solution (3) steady state response for 
longer baselines. 

Since the geometric solutions are not affected by the erroneous 
modeling of either the satellite motion or the motion of reference frames, 
we have taken solution (3) to constitute the standards of comparison in 
assessing the effectiveness of the SRD method versus the dynamic mode 
methods. Solution (3) was chosen over the corresponding minimum 
constrained solutions because, of all of the available events, the steady 
state response for more than the two very short baselines (i.e.. 
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7110-7220 and 7109-7886) was not possible via the minimum constrained 
solutions (Table 8). 

The only errors affecting solution (3) result from the ineffectiveness 
to implement the ill-defined scale by fixing the third coordinate of 
station 7122 to that of the (CSR)85L01 solution (see Section 4.4.2). This, 
however, is difficult to assess, but on the basis of the results presented 
in the next section it can be safely stated that the scale along a band 
in the direction of baseline 7109-7122 has been properly defined by 
constraining this baseline. In fact the value of the scale factor varies 
from solution C2 to that of C3, from 6.61x10“' to 7.72x10“' for baselines 
longer than 627 km while for the two baselines of about 274 km the 
value of the scale factor drops down to 4.5x10“'. This dependency of 
the scale factor on the baseline lengths is expected since for shorter 
baselines the implied geometry is stronger. 


4.4 BASELINE ESTIMATION VIA THE STEADY STATE RESPONSE OF THE 
SRD METHOD 

The incorporation of additional observations into the SRD solutions will 
lead to a steady state response for the same reasons mentioned in the 
previous section. This response is even faster if observations of 
"improved quality" are available. The term improved quality indicates 
improvement of the weight coefficient matrix. For instance, from two 
sets of observations (a) and (b) having identical sizes, the set (b) is of 
improved quality if the difference (P^ - P,) of their weight matrices is 
a positive (semi) definite matrix. Based on this, it is a trivial exercise to 
prove that the matrix (Qx* - Qx**) is positive (semi) definite if the design 
matrix is the same for both sets of observations (a) and (b) and of full 
rank. Qjj* and QJ** denote the weight coefficient matrices of the same 
adjusted parameter vector obtained on the basis of sets (a) and (b). 
Since (Q$® - Qx**) is a positive (semi) definite matrix, it follows that 

TrCQ^b) < TrCQja) (4-2) 

Although the two sets (a) and (b) are drawn from different populations, 
consistency between the accuracy of the models employed and the 
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accuracy of the observations leads to a close-to-unity a posteriori 
variance of unit weight, thereby reducing, as it follows from equation 
(4-2), the a posteriori variances of the adjusted parameter vector. This 
vector contains the earth-fixed Cartesian coordinates of the ground 
stations and the inertial initial state vectors which are treated in the 
present study as nuisance parameters. The adjusted station coordinates 
and their statistics are used to estimate the baseline lengths and their a 
posteriori standard deviations. Reduced variances in station 

coordinates, obtained from improved quality observations, yield baselines 
the variances of which are also reduced only if the non-linearity of the 
employed models allows such a reduction to take place. These 
reductions, taking place for most of the SRD solutions presented in the 
next section, imply a faster steady state response. 

Constraining of estimable quantities in any SRD solution results in 
baselines with reduced variances thereby implying again a faster steady 
state response. Thus, incorporation of additional observations, improved 
quality observations, and constraining of estimable quantities constitute 
the three major factors leading to a steady state response. In the SRD 
solutions, this response is achieved by balancing the contribution of the 
first and the third factors because we have no control over the second 
factor since the quality of already recorded observations does not 
change. Furthermore, the goal has been to restrict as much as possible 
the contribution of the third factor because the steady state response 
achieved on the basis of constraining estimable quantities may be 
affected severely by the errors affecting those quantities. These 
quantities are in error because their recovery has also been based on 
erroneous observations. If the propagation of the errors affecting the 
constrained estimable quantities is not well controlled, it will lead to 
erroneous recovery of the adjusted parameter vector and therefore to 
erroneous baseline estimates. Even worse, the a posteriori variances of 
those estimates might also be reduced making it more difficult to 
reliably assess their accuracy. 
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4.4.1 Pass Baseline Geometry and Its Manifestation in the Design Matrix 
of the SRD Observable 

The underlying characteristic inherently present in any factor leading 
to a steady state response is its ability to strengthen the geometry 
implied in the SRD solution under question. Thus, with stronger 
geometry, the steady state response will be achieved on the basis of a 
reduced number of observations. In other words, observations 

containing more information for the recovered estimable quantities will 
lead to a faster steady state response. 

Pig. 17 shows the horizontal plane of the starting point (1) of the 
baseline (12). The ending point (2) of the baseline and the 

simultaneously observed satellite positions have been projected on this 
plane. This simplification is aimed at revealing the information the SRD 
observables contain about the estimated baseline for two characteristic 
geometric configurations, namely, one when the subsatellite tracks are 

parallel to the estimated baseline (Pig. 17a) and one when the 

subsatellite tracks are perpendicular to that baseline (Pig. 17b and c). 

It is a matter of trivial trigonometric manipulations to confirm that 
with subsatellite tracks parallel to the estimated baseline the magnitude 
of the SRD observable tends towards the length of the estimated 
baseline as the subsatellite point moves away from either end of the 
baseline along the subsatellite track denoted by (e) in Pig. 17a. Thus, 
with this geometry the SRD observable directly relates to the length of 
the baseline. However, as the subsatellite point moves towards the 
perpendicular bisector of baseline (12'), the SRD observable tends 
towards zero. 

With the subsatellite tracks perpendicular to the estimated baseline, 
the SRD observable tends towards zero as the subsatellite point moves 
away from the projected baseline along the line (c) in Pig. 17b and c. 
This observable also tends towards zero as the subsatellite track (e) 
tends toward the perpendicular bisector of the projected baseline (12'). 
However, the SRD observable tends towards the baseline length as the 
subsatellite point moves toward point (I) and as this point moves in the 
direction of the baseline and away from it (Pig. 17c). Thus, for a faster 
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steady state response, parallel passes, the coobserved parts of which 
are extended well away from the endpoints of the baseline, should be 
preferred. When perpendicular passes are available, the ones located 
well away from the perpendicular bisector plane of the estimated 
baseline should be preferred. 

The previously described favorable geometry manifests itself in the 
strong independence prevailing among the columns of the resulting 
design matrix. The term "strong independence" indicates the sharp 
variation characterizing the ratio of the corresponding entries between 
any two columns of the design matrix. The sharp variation resulting 
from favorable geometry is seen in the partial derivatives of the SRD 
observable taken with respect to the coordinates of one of the ground 
stations and to those of the inertial initial state vector (Appendix A) 

(4-3) 

Pj2 



Xj - Y2 

Pj2 


Xj - Y, 


T 

(SNP) 



(4-4) 


where Xj, Y,, Y 2 are the earth-fixed position vectors of the satellite at 
epoch j and stations (1) and (2) respectively. 

Since the partial derivatives of equation (4-3) have a common 
denominator, the variation of their ratios is controlled by the variations 
prevailing among the coordinates of the satellite as it moves along the 
coobserved part of the pass referred to from now on as "common part." 
The longer the common parts, the sharper the variations are for the 
ratios of the partial derivatives of equation (4-3). Common parts tend 
to be longer as the baseline decreases thereby resulting in sharper 
variations among the ratio of those partial derivatives for shorter 
baselines. However, these variations become sharper for both short and 
long baselines as more and more passes of different spatial distribution 
are incorporated into the solution. For long baselines (i.e., > 3500 km), 
the common part tends to be short if the pass is parallel to the baseline 
thereby reducing the variation of the ratio of those partial derivatives. 
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For perpendicular passes, the common parts tend to be a little longer 
for long baselines, allowing therefore for sharper variations in the ratio 
of the partial derivatives of equation (4-3). 

The variations of the ratios between the partial derivatives of 
equation (4-4) as well as between the partial derivatives of equations 
(4-3) and (4-4) tend for perpendicular passes to be substantially 
reduced. The variations of the ratios among the partial derivatives of 
equation (4-4) are controlled by the changes in the observed ranges pja 
and Pji, the changes in the satellite coordinates and the changes in the 
columns of the state transition matrix. When the SRD observables tend 
to become shorter (Fig. 17b and c), the observed ranges tend to be 
equal and therefore the coordinates of the coobserved satellite positions 
tend to cancel out. With such cancellations taking place and with the 
tendency to have equal denominators in those partial derivatives, the 
variations of their ratios are primarily controlled by those among the 
transformed state transition matrix ((SNP)’3Rj/3Xo, ...), the columns of 
which, apart from the common denominator, tend in this case to be 
multiplied by the components of the estimated baseline. Thus, for 
shorter baselines the variations among these partial derivatives tend to 
be sharper than those for longer baselines (Appendix B). The 
unfavorable geometry associated with almost perpendicular passes 
manifests itself with an almost constant numerator in the partial 
derivative taken with respect to the station coordinate measured along 
the direction of the estimated baseline. This happens because in this 
direction the satellite coordinates change very little as the satellite 
moves along a pass almost perpendicular to the direction of the baseline. 
This geometry also leads to small variations of the numerators in the 
partial derivatives of equation (4-4); therefore, in the design matrix the 
ratio between the entries of the columns corresponding to the station 
coordinate, recorded along the baseline direction, and to the coordinates 
of the initial state vectors will exhibit very little variation because the 
denominators of the corresponding entries are nearly equal. 

Reduced relative variations among the columns of the design matrix 
lead to high correlations among the corresponding adjusted parameters. 
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These correlations tend to increase or decrease as the homogeneity and 
the strength of the implied geometry decreases or increases 

respectively, a fact also manifested in the resulting poor or strong 
conditioning of the normal matrix. 

The above discussion, the aim of which is to help in understanding 
the results presented in the next section, evidently reveals that in 
single SRD baseline solutions, the strength of the geometry fades away 
as the length of the estimated baseline increases. For longer baselines 
(i.e., > 4000 km), the steady state resi>onse cannot be reached because 
poor geometry results in an ill-conditioned normal equation matrix. 

The weak geometry associated with long single baseline solutions is 
substantially improved through a network solution. Such a solution has 
not been performed in the present study because of the large number 
of observations and the limiting computer capabilities. If SRD normal 
points were available a network solution would have been possible 
which, however, is beyond the scope of this study. 


4.4.2 Baseline Results 

This section describes the process of achieving steady state response 
for those baselines only for which such a response, through single SRD 
baseline solutions, was possible on the basis of the data collected during 
the main MERIT campaign. 

The adjusted parameter vector in single SRD baseline solutions 


contains the coordinates of the baseline end-points and the components 
of the initial state vectors for all of the arcs involved. In these 
solutions, obtained via a weighted least squares adjustment, the 
coordinates of one baseline end are held fixed while the coordinates of 
the other end and the components of the initial state vectors are 
allowed to adjust by assigning, through their weight matrix, the 
following standard deviations 

= 0.0001 m (fixed baseline end) 

("free” baseline end) 

= = 50 m (initial position) 

(initial velocity) 




<Tx ffy 


= <T, 


<Tj; = 20 m 




= 5 cm/s 
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The term "free" designates the adjusting baseline end. 

The standard deviation of 20 m assigned on the coordinates of the 
free baseline end reflects an accuracy estimate of their approximate 
values. These values have been obtained by altering from 20 m to 
100 m the coordinates of the (CSR)85L01 solution (Tapley et al., 1985b). 
The standard deviations of 50 m and 5 cm/s assigned to the approximate 
values of the components of the initial state vectors, reflect the 
accuracy of those values. They have been obtained by integrating, over 
a period of a month, Lageos' equations of motion using as initial state 
vectors those computed at EG&G for the entire main MERIT campaign 
(Pavlis, 1986; private communication). At the beginning of each month 
the EG&G estimates differ from those obtained by integrating Lageos’ 
equation of motion over the period of the previous month because a 
simplified orbital modeling was used in the integration. These 
differences ranging from 50 m to 150 m for the position and from 5 cm/s 
to 13 cm/s for the velocity reflect the assigned accuracies to the 
components of the initial state vector. 

Constraining one baseline end, by applying large weights to its 
coordinates, results in implicit constraining of its latitude and geocentric 
distance which both constitute, in the dynamic/semidynamic environment, 
a set of estimable quantities if the scale has been defined through the 
appropriate adoption of units and constants (i.e., units of time, velocity 
of light, etc.) and if a solution is possible through appropriately 
adopting physical constraints (Van Gelder, 1978). Weighting moderately 
the components of the initial state vectors results also in moderate 
implicit weighting of the estimable quantities associated with the 
geometric characteristics of the orbit. Adjusting only one baseline end 
and the initial state vectors, in either dynamic or semidynamic single 
baseline solutions, results in implicit and/or explicit constraining of a 
large number of estimable quantities. For instance constraining the 
and potential coefficients together with the geocentric gravitational 

constant (GM^) and the mean radius of the earth ( gl ^) results in explicit 
constraining of the A„^ and coefficients which in the dynamic/ 

semidynamic environment constitute a set of estimable quantities (ibid.) 
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Erroneous constraining of estimable quantities distorts the geometric 
and/or physical characteristics of the dynamic environment, thus if 
enough time is allowed the accumulated errors can exceed the noise level 
of the observations to the extent that the adjusted estimable quantities 
will be affected by errors larger than those implied by the accuracy of 
the observations. Even worse with a priori variance of unit weight 
equal to one, the a posteriori standard deviations of these quantities 
will also be reduced only if the a posteriori variance of unit weight is 
close enough to unity so it preserves the relative structure of the 
weight coefficient matrices (ibid.) 

Under those conditions baselines estimated on the basis of Cartesian 
coordinates will have their a posteriori variances reduced only if the 
non-linearity of the models employed allow such reductions to take 
place. In single 3RD solutions the errors affecting the baselines, 
computed from the erroneously adjusted Cartesian coordinates of the 
free baseline end, result from the errors originated by erroneously 
constraining the parameters entering in the computation of the 
discrepancy vector L of equation (2-60). Erroneous baseline lengths 
with reduced variances constitute a potential source for misleading 
inferences in regard to the accuracy of the estimated baselines. 
Therefore, the a posteriori standard deviations of the estimated 
baselines will only be used , as in the geometric solutions, to indicate 
how far away the solution is from its steady state. The actual accuracy 
of the estimated baselines will be inferred on the basis of the 
comparison with the baselines estimated through geometric method 
(present study) and range dynamic methods (CSR at UTX and ZIPE at 
Potsdam) (next section). 

Baselines estimated through Cartesian coordinates will also be 
Eiffected in the dynamic environment by the errors commited in the 
implementation of the terrestrial reference frame, as it happens, for 
instance, when the barycenters of the observations of the baseline end 
stations are well apart in time (Pavlis et al., 1983). Strict simultaneity, 
however, implied by the nature of the 3RD observable eliminates such 
inconsistencies in the present study. 
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For various baseline lengths ranging from 8 m to 4000 km, Tables 9 
through 17 list the final results and those obtained at intermediate 
stages of the process leading to a steady state response. The a priori 
variance of unit weight is taken equal to one for all of the baseline 
solutions presented in this section. Fig. 18 shows the locations of the 
stations involved in the SRD solutions (lower part) together with a 
typical LAGEOS groundtrack for some of their coobserved passes (upper 
part). 

The lengths of these common parts increase or decrease as the 
length of the baseline being estimated decreases or increases 
respectively. 

Tables 9 and 10 contain the results of the steady state response 
reached for two very short baselines via both short and long arc 
modes. The first column lists the number of passes coobserved by the 
baseline end points while the second lists the number and the duration 
of the arcs, the position and orientation of which were adjusted to fit 
"best" the available observations. For instance, 6 (Ih) in the first row 
of Table 9 indicates that six arcs were adjusted, each of one hour long, 
and 1 (4h) and 1 (2d) in the fifth row indicate that observations from 
eight passes (first column of this row) were adjusted using two arcs, 
one four hours long, and one two days long. The third, fourth and 
fifth column list the total number of the SRD observables, the a 
posteriori variance of unit weight and the root mean square of all of the 
SRD residuals respectively. The last column contains the estimated 
baseline length together with its a posteriori standard deviation. All of 
the remaining tables through 17 have the same format except for some 
obvious very minor differences. 

For baselines 7109 - 7886 and 7110 - 7220 steady state response has 
been reached via both short arc and long arc solutions on the basis of 
10 passes (Table 9, rows 3-7) and 17 passes (Table 10, rows 2-5) 
respectively. As in the geometric solutions, the steady state response is 
also associated with a posteriori standard deviations below the 
centimeter level. After the steady state, the estimated lengths of these 
two baselines do not change at the centimeter level either through a 
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Fig. 18 Station location and Lageos groundtracks 
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Table 9 Steady State Response of Baseline 7109 - 7886, Parallel 
Passes 

All length units in meters. 


No. of 
Passes 

Integration 
Lengths ^ * ) 

No. of 
Observ. 


RMS of All 
Residuals 

SRD Baseline 

7. 

Short— Arc McxIg 







6 

6 (1 h) 

8,414 

0.94 

0.070 

.719 

± 

0.002 

8 

8 (1 h) 

12,807 

0.90 

0.069 

.731 

± 

0.002 

10 

10 (1 h) 

18,589 

0.86 

0.067 

.737 

± 

0.001 

13 

13 (1 h) 

25,865 

0.85 

0.067 

.738 


0.001 

Long-Arc Mode 







8 

1 (4 h) 

13,527 

0.90 

0.069 

.733 

± 

0.002 


1 (2 d) 







10 

1 (4 h) 

20,216 

0,86 

0.067 

.741 

± 

0.001 


1 (4 d) 







13 

1 (4 h) 

27,697 

0.85 

0.067 

.738 

± 

0.001 


1 (7 d) 







(0 

k (i h) = k 

arcs of i 

hours 






k d) = k 

arcs of 1 

days 





(O 

a posteriori 

variance 

of unit 

weight 





short arc or a long arc solution and furthermore these lengths are the 
same for both short arc and long arc solutions, thereby indicating that 
the accumulated orbital biases for so short baselines cancel out 
completely, a plausible property of the SRD observable. The correlations 
among the components of the adjusted parameter vector are for so short 
baselines substantially reduced because the lengths of the coobserved 
parts of the passes tend to be long and their orientation is not 
important since the subsatellite tracks for these short baselines are 
located well away from their end points (Section 4.4.1). 

The 34 cm a posteriori standard deviation of the baseline solution, 
listed in the first row of Table 11, shows that although 17 passes of 
about one hour long are available, steady state response, with a short 
arc solution, is not possible for baseline 7110 - 7265 because the 
pass-baseline geometry has deteriorated as the length of the estimated 
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Table 10 Steady State Response of Baseline 7110 - 7220, Parallel 
Passes 

All length units in meters. 


No. of 
Passes 

Integration 
Lengths ( ^ ) 

No. of 
Observ. 

( 2 ) 

<T 0 

RMS of All 
Residuals 

SRD Baseline 
15. 

Sbort-Arc Mode 








11 

11 

(1 

h) 

12,040 

0.95 

0.084 

.233 

± 

0.003 

17 

17 

(1 

h) 

17,971 

1.09 

0.091 

.236 

± 

0.002 

Long-Arc Mode 








17 

1 

(4 

h) 

17,982 

1.07 

0.090 

.240 

± 

0.002 


1 

(1 

d) 








1 

(2 

d) 








3 

(3 

d) 







17 

1 

(1 

h) 

17,983 

1.08 

0.090 

.250 

± 

0.002 


1 

(3 

d) 








1 

(4 

d) 








2 

(5 

d) 







17 

1 

(1 

h) 

17,984 

1.08 

0.090 

.238 

± 

0.002 


1 

(3 

d) 








1 

(6 

d) 








1 

(7 

d) 








(i) k(ih)=k arcs of i hours 
k (i d) = k arcs of i days 


( 2 ) a posteriori variance of unit weight 

baseline has increased from 8 m (Tables 9) to about 274 km (Table ll). 
This deterioration is even worse since the geometry of the orbit is such 
that passes parallel to this baseline do not exist and only passes 
intersecting it at about * 30 to * 50 degrees are available creating, 
therefore, a geometry which is worse than that of the parallel passes 
and better than that of the perpendicular ones. Steady state response, 
however, for this solution can be reached by strengthening its 
geometric characteristics on the basis of constraints imposed on some 
additional estimable quantities. In the dynamic environment this is 
accomplished by increasing the maximum length of continuous integration 
(maximum arc length) which intensifies, through the implied geometric 
strength, the effect of the coefficients An„, and Bn„, on the resulting 
long orbital arcs. Care should be excercised not to destroy the 
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Table 11 Steady State Response of Baseline 7110 - 7265, Passes 
Within *30* - ±50* 

All length units in meters. 


No. of Integration No. of a*) RMS of All SRD Baseline 

Passes Lengths Observ. Residuals 2740 . 


Sbort-Arc Mode 


17 

17 (1 h) 

21,767 

0.77 

0.0734 

70.453 

± 

0.335 

Long-Arc Mode 







17 

1 (1 d) 
3 (2 d) 

21,780 

0.80 

0.075 

69.391 

± 

0.012 

17 

2 (2 d) 
1 (4 d) 

21,781 

0.80 

0.075 

69.482 

± 

0.009 

17 

1 (2 d) 

21,781 

0.80 

0.075 

69.494 

± 

0.009 


1 (3 d) 
1 (7 d) 


(1) k(^h)=k arcs of £ hours 
k (i d) = k arcs of t days 

( 2 ) a posteriori variance of unit weight 

accuracy of the available observations by the errors accumulated, over 
these long arcs, due to errors affecting the and B„n, potential 

coefficients. 

With this in mind and on the basis of all observations, three long 
arc solutions have been performed by allowing maximum arc lengths of 
up to three, five and seven days respectively. The results of these 
three solutions are listed in 2nd through 4th row of Table 11. In the 
first long arc solution, employing one arc of one day long and three 
arcs of two days long, the length of the estimated baseline changed by 
106 cm as compared to that of the short arc solution. Steady state 
response, however, for this solution has not yet been reached since the 
length of this baseline changes by about 10 cm when arcs up to five 
and seven days are allowed (Table 11, rows 3-4). However, increasing 
the maximum arc length from four to seven days (rows 3-4) results in a 
change of the estimated baseline length of about 1 cm, thereby implying 
that steady state response is being reached and that the baseline length 
of the solution listed in 3rd row of this table is the least affected by 
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Table 12 Steady State Response of Baseline 7109 - 7110, Parallel 
Passes 

All length units in meters. 

No. of Integration No. of 2 ( 2 ) RMS of All SRD Baseline 

Passes Lengths Observ. Residuals 883602. 


Sbort-Arc Mode 


4 

4 

(1 

h) 

11,512 

0.5 

0.029 

.676 

* 0.30 

8 

8 

(1 

h) 

36,236 

0.36 

0.026 

.309 

* 0.04 

12 

12 

(1 

h) 

59,011 

0.35 

0.026 

.249 

* 0.03 

16 

16 

(1 

h) 

69,083 

0.43 

0.026 

.251 

* 0.03 

Long-Arc Mode 







12 

1 

(1 

h) 

59,020 

0.36 

0.026 

.220 ± 

0.001 


1 

(2 

d) 







1 

(3 

d) 






12 

1 

(4 

h) 

59,021 

0.36 

0.026 

.224 * 

0.001 


1 

(3 

d) 






12 

1 

(7 

d) 

59,022 

0.86 

0.042 

.217 * 

0.002 

16 

2 

(2 

d) 

69,816 

0.37 

0.027 

.225 ± 

0.001 


1 

(3 

d) 






16 

1 

(2 

d) 

69,817 

0.38 

0.027 

.226 * 

0.001 


1 

(4 

d) 






16 

1 

(1 

d) 

69,817 

0.79 

0.057 

.219 * 

0.001 


1 

(6 

d) 







(1) k(^h)=k arcs of i hours 
k d) = k arcs of t days 

( 2 ) a posteriori variance of unit weight 


any accumulated orbital errors, because at steady state this solution 
employs the shorter long arcs and the a posteriori variance of unit 
weight is smaller than one. The correlations among the components of 
the adjusted parameter vector do not cause, for this baseline, any 
instability since almost all of them are less than 0.90 with very few just 
exceeding this value. 

Tables 12 and 13 contain the results for the steady state response 
of two baselines with the same length, the same orientation, and one 
common end point occupied by station 7110. The observations recorded 
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Table 13 Steady State Response of Baseline 7110 - 7886, Parallel 
and Perpendicular Passes 
All length units in meters. 


No. of Integration 
Passes Lengths ( ‘ ) 

No. of 
Observ, 

( 2 ) 

<r 0 

RMS of All 
Residuals 

SRD Baseline 
883606. 

Sbort-Arc Mode 






33 33 (1 h) 

parallel 
passes 

58,261 

0.74 

0.064 

.467 * 

0.056 

33 39 (1 h) 

parallel 
+ 6 

perpendi cul ar 
Long-Arc Mode 

61,037 

0.75 

0.064 

.459 ± 

0.056 

33 4 (1 h) 

parallel 4 (4 h) 

passes 2 (1 d) 

1 (1.5 d) 
3 (2 d) 

58,248 

0.75 

0.065 

.347 ± 

0.004 

33 2 (1 h) 

parallel 1 (4 h) 

passes 1 (Id) 

3 (2 d) 
3 (3 d) 

58,252 

0.76 

0,065 

.342 * 

0.003 

33 1 (4 h) 

parallel 1 (2 d) 

passes 1 (3d) 

1 (4 d) 

1 (4.5 d) 
1 (7 d) 

58,256 

0.77 

0.065 

.335 i 

0.003 


(i) k(^h)=k arcs of i hours 

k (i d) = k arcs of t days 


( 2 ) a posteriori variance of unit weight 

by stations 7109 and 7886, occupying the other ends of these two 
baselines, have an accuracy of 0.028 m and 0.070 m respectively, 
therefore, examination of the results presented in these two tables will 
show how the quality of the observations affects the speed of the 
steady state response (Section 4.4). Steady state response of baseline 
7109 - 7110 has been reached via short arc solution on the basis of 12 
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passes having approximately 59,000 observations, as is confirmed by 
comparing the third and fourth rows of Table 12. The a posteriori 
standard deviation of the recovered baseline is at the 3 cm level 
because, in short arc solutions, the components of the initial state 
vectors are recovered with a relatively low precision (20 m to 50 m). 

To examine the effects of using long arc solutions to achieve steady 
state response, we have performed six long arc solutions, three of 
which are based on observations from 12 passes (Table 12, rows 5-7) 
and three on observations from 16 passes (rows 8—10) when maximum arc 
lengths up to three, five and seven days are allowed respectively. The 
baseline lengths, estimated through all of these long arc solutions, do 
not change their length at the centimeter level, when either the number 
of observations or the maximum arc length increases and furthermore, 
the a posteriori standard deviation of the estimated baselines assumes 
millimeter level values. However, the baseline length changed from the 
short arc solution to the long arc solution by about 2.5 cm due to 
accumulation of orbital errors. The RMS of all of the SRD residuals for 
the long arc solutions (Table 12, rows 5-6,and 8-9) is about the same as 
that of the short arc ones (rows 1-4), but it is twice as much when arcs 
up to seven and six days are allowed in those solutions (rows 7 and 10). 

Steady state reached through a short arc solution assumes a 
posteriori standard deviation of 3 cm (Table 12), thus, for baseline 
7110-7886 shown in Table 13 such a response has not yet been reached, 
because a short arc solution, on the basis of all of the available 
observations, results in a 6 cm a posteriori standard deviation of the 
estimated baseline length (Table 13, row 1). This claim is also confirmed 
by the 10 cm change of the estimated baseline length when long arc 
solutions are performed (rows 3-5). The perpendicular passes 
incorporated in the short arc solution (2nd row), tend, in the light of 
weak geometry, to decrease the estimated baseline length because for 
those passes the SRD observable tends to become very short (Section 
4.4.1). Thus, the 8 mm decrease (2nd row) on the basis of just 2,776 
additional observations is another indicator confirming that steady state 
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response for this baseline through a short arc solution has not been 
reached. Consequently, such a response is sought through long arc 
solutions performed on the basis of arcs up to two, three, and seven 
days long (Table 13, rows 3-5). Since in these three long arc solutions 
the estimated baseline length did not change at the cm level and since 
its a posteriori standard deviation is below the cm level, it is assumed 
that steady state has been reached in those solutions. 

For this baseline, although of same length and direction as baseline 
7109 - 7110, steady state response was not possible on the basis of 33 
passes with 58,261 observations processed via a short arc solution, 
because the accuracy of the observations of station 7886 is twice as 
large as those of station 7109. This confirms the claim, made in Section 
4.4, that steady state response on the basis of improved quality 
observations is faster. 

As of the correlations among the components of the adjusted 

parameter vector, they are well below the 0.90 value for all short arc 
solutions shown in Tables 12 and 13, except for very few of them, 
among the components of some of the recovered state vectors, that 
tend to be Just a little higher than this value without, nevertheless, 
affecting the conditioning of the normals to the extent that could result 
in an algorithmically singular normal equations matrix. This behavior of 
the correlations is not surprising since in estimating these two baselines 
only parallel passes, the coobserved parts of which are extended well 
beyond the baseline end points, were employed, thereby, leading 

according to the discussion in Section 4.4.1 to a favorable geometry 
which is manifested in the reduced correlations among the components of 
the adjusted parameter vector. On the contrary, when in some of the 
long arc solutions shown in Tables 12 and 13, arcs of one or four hours 

long are employed together with arcs of two, three or seven days long, 

the relative geometry of the hours long arcs is very weak as compared 
to that implied by the days long arcs, thereby leading to high 
correlations among the components of initial state vectors of the hours 
long arcs. This inhomogeneity together with the potential to have many 
passes with weak geometry, when long baselines of about 4000 km are 
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estimated^ results in an algorithmically singular normal equations matrix 
(see below). 

Baseline 7110 - 7086, the steady state response of which is shown in 
Table 14, has a relatively large number of coobserved passes each 
having a substantially reduced number of observations as opposed to 
those recorded in the passes used in the estimation of the baselines 
described in Tables 9 through 13. Examination of the steady state 
response for this baseline will reveal the effect of the geometry, implied 
by the many passes, as opposed to the number of the available 
observations. The 41 cm a posteriori standard deviation obtained with 
the short arc solution, shown in the first row of this table, indicates 
that steady state has not yet been reached. Thus, six long arc 
solutions have been performed, the three of which use observations from 
29 passes (Table 14, rows 2-4) while the other three use observations 
from 40 passes (rows 5-7) and arc lengths of up to three, five and 
seven days respectively. On the basis of 29 passes, processed in the 
long arc mode, steady state response has not been reached because 
changing the maximum arc length from five days to seven days the 
length of the estimated baseline changes as much as 18 cm (Table 14, 
rows 3-4). However, when 11 more passes with about 4,000 additional 
observations are included, the lengths of the estimated baseline change 
only at the centimeter level (rows 4-7) and their a posteriori standard 
deviations has dropped below the centimeter level, thereby suggesting 
that steady state response for this baseline has been reached. Since 
the a posteriori variance of unit weight is smaller than one, the 
estimated lengths are assumed not to have been influenced by any 
errors accumulated over these long periods and therefore, the length 
estimated with the smaller standard deviation is assumed to be the 
closest to that of the steady state (last row of this table). 

In the short arc solution (1st row of Table 14), the correlations 
among the components of the adjusted parameter vector are almost all 
of them less than 0.90 with few just exceeding this value. These 
correlations (i.e., >0.90) exist among some components of the initial state 
vectors corresponding to arcs, the geometry of which is not favorable 
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Table 14 Steady State Response of Baseline 7110 - 7086, Passes 
Within *20* - ±60* 

All length units in meters. 


No. of 
Passes 

Integration 
Lengths ( * ) 

No. of 
Observ . 

<^o 

RMS of All 
Residuals 

SRD Baseline 
119829_. 

Sbort-Arc Mode 






40 

40 (1 h) 

16,752 

0.76 

0.078 

1.110 ± 

0.411 

Long-Arc Mode 






29 

8(lh), l(4h) 
3(ld), 2(2d) 
3(3d) 

12,395 

0.81 

0.081 

0.814 ± 

0.042 

29 

5(lh), 2(4h) 
2(ld), 2(2d) 
2(4d), 2(5d) 

12,414 

0.82 

0.081 

0.821 i 

0.036 

29 

4(lh), l(ld) 
l(2d), 2(4d) 
2(5d), 2(7d) 

12,400 

0.84 

0.083 

0.998 ± 

0.015 

40 

lO(lh), l(4h) 
3(ld), 2(2d) 
5(3d) 

16,771 

0.82 

0.081 

0.987 ‘ 

0.009 

40 

7(lh), l(4h) 
3(ld), 3(2d) 
l(3d), 2(4d) 
2(5d) 

16,790 

0.80 

0.080 

1.024 i 

0.009 

40 

7(lh), l(ld) 
l(2d), l(4d) 
2(5d), 2(6d) 
2(7d) 

16,776 

0.87 

0.084 

1.005 ± 

0.008 


(i) k(ih)=k arcs of i hours 
k (-^ d) = k arcs of £ days 


( 2 ) a posteriori variance of unit weight 

with respect to the estimated baseline (Section 4.4.1), that is, when the 
passes are close to being perpendicular to the baseline and/or when 
their coobserved part is short and close the perpendicular bisector 
plane of the estimated baseline. 

Increasing the maximum arc length results in arcs of different 
lengths ranging from one hour to either three, five, or seven days long, 
all of which are employed in the same long arc solution. This in turn 
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weakens the geometric strength of the short arcs because in those 
solutions the long arcs are the ones dominating the implied geometric 
strength. The weaker the geometry of those short arcs, in a short arc 
solution, the higher the correlations among their components adjusted in 
long arc solutions, in which these short arcs could not be matched with 
any other arcs. Furthermore, if only passes of weak geometry are 
matched together they result in a long arc also of weak geometry, 
thereby leading again to high correlations among the components of its 
adjusted initial state vector. These correlations become larger as the 
maximum arc length increases and the long arcs of weak geometry 
cannot be matched with any other ones of strong geometry. This 
pattern, which is present in the long arc solutions of Table 14, does 

not lead to an algorithmically singular normal equation matrix, simply 
because the length of the baseline is not long enough to result in a 
geometry so weak that could lead to near singularities, although passes 
parallel to this baseline do not exist (Fig. 18). 

Table 15 contains the results for the steady state response of 
baseline 7110 — 7122, the geometry of which is such that passes 
parallel to this baseline do exist and in fact all of the passes employed 
in the solutions of Table 15 are parallel to this baseline. Since steady 
state response for this baseline was not possible through a short arc 
solution, long arc solutions have been performed on the basis of 15 
passes (rows 1-2) and 22 passes (rows 3—5) with maximum arc lengths 
up to three, five and and seven days respectively. The temporal 
distribution of the observed passes is such that when 15 of them are 
used the solution allowing maximum length of three days coincides with 
that allowing maximum arc length of five days, and therefore only one of 
them is shown in Table 15 (1st row). These passes contain observations 
from station 7122 before and after it was upgraded (Table 2). 

The 3 mm and 2 mm level of the a posteriori standard deviations 
associated with the baselines obtained through all of the long arc 
solutions, shown in Table 15, indicate that a steady state response for 
this baseline has been reached. This is also confirmed by the fact that 
the estimated baseline length changes at just the centimeter level when 


155 



Table 15 Steady State Response of Baseline 7110 - 7122, Long-Arc 
Mode, Parallel Passes (all length units in meters) 


No. of 
Passes 

Integration 
Lengths ( ‘ ) 

No. of 
Observ . 

0 

RMS of All 
Residuals 

SRD Baseline 
1437139. _ 

15 

1 (1 h) 
5 (1 d) 
1 (2 d) 

42,328 

0.65 

0.044 

.307 ± 0,003 

15 

1 (1 h) 
4 (1 d) 
1 (6 d) 

42,329 

0.66 

0.045 

.309 * 0.002 

22 

2 (1 h) 
4 (1 d) 

3 (2 d) 
1 (3 d) 

69,803 

0.52 

0.035 

.305 * 0.002 

22 

1 (1 h) 
5 (1 d) 

2 (2 d) 
1 (3 d) 

69,804 

0.54 

0.037 

.302 * 0.002 

22 

2 (1 d) 
1 (2 d) 

1 (6 d) 

2 (7 d) 

69,807 

0.81 

0.055 

.293 t 0.002 

(O 

k (i h) = k 
k (i d) = k 

arcs of i 
arcs of i 

hours 

days 



( 2 ) 

a posteriori 

variance 

of unit 

weight 


on the 

basis of 22 passes the 

maximum arc lengths 

change from 3 


days (Table 15, rows 4-5), The correlations follow the same pattern as 
that described for baseline 7110 - 7086 but since the passes are now 
parallel to the estimated baseline, weak geometry is implied only by 
those the coobserved part of which is extended only in between the 
baseline endpoints (Section 4.4.1). However, the geometry implied by 
those passes is not so weak to cause any algorithmic singularity in the 
normal equations matrix. Such weakness in the geometry causes 
problems when the lengths of the estimated baselines increases to 
3500 km and 3700 km (Tables 16 and 17). 
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Table 16 Steady State Response of Baseline 7109 - 7105, Long- Arc 
Mode, Passes Within *20* - *30* 

All length units in meters. 


No. of Integration 
Passes Lengths ( * ) 

No. of 
Observ. 

<T0 

RMS of All 
Residuals 

SRD Baseline 
3703351. 

Maximum Arc Length = 

2 days 




42 

10 (Ih), 5(4h) 
2(7h), 6(ld) 

78,801 

0.54 

0.032 

.628 * 0.003 

54 

12(lh), 5(4h) 
2(7h), 8(ld) 

106,770 

0.52 

0.032 

.686 * 0.003 

62 

14(lh), 5(4h) 
2(7h), 9(ld) 

129, 147 

0.50 

0.031 

.702 * 0.002 

72 

18(lh), 5(4h) 
2(7h),ll(ld) 

151,901 

0.49 

0.031 

.693 * 0.002 

Maximum Arc Length = 

3 days 




42 

5(lh), 2(4h) 
l(7h), 4(ld) 
2(2d), 3(3d) 

72,783 

0.57 

0.033 

.677 * 0.003 

54 

5(lh), 2(4h) 
l(7h), 4(ld) 
4(2d), 4(3d) 

106,989 

0.58 

0.033 

.689 * 0.002 

62 

5(lh), 2(4h) 
l(7h), 5(ld) 
5(2d), 4(3d) 

129,155 

0.54 

0.032 

.699 * 0.002 

72 

6(lh), 2(4h) 

151,910 

0.53 

0.032 

.695 * 0.002 


l(7h), 7(ld) 
6(2d), 5(3d) 
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Table 16 (cont’d) 


No. of 
Passes 

Integration 
Lengths ( * ) 

No. of 
Observ. 

( 2 ) 

/r2 

^ 0 

RMS of All 
Residuals 

SRD Baseline 
3703351. 

Maximum Arc Length = 

7 days 




42 

4(lh), l(4h) 
l(7h), l(ld) 
l(2d), l(3d) 
l(5d), 4(6d) 
l(7d) 

78,809 

1.11 

0.046 

.706 * 0.004 

54 

4(lh), l(4h) 
l(7h), l(ld) 
l(2d), l(3d) 
l(5d), 5(6d) 
2(7d) 

106,780 

3.33 

0.081 

.815 * 0.005 

62 

4(lh), l(4h) 
l(7h), l(ld) 
l(2d), l(3d) 
l(5d), 6(6d) 
2(7d) 

129,159 

2.78 

0.073 

.809 * 0.005 

72 

4(lh), l(4h) 
l(7h), 2(ld) 
l(2d), l(3d) 
l(4d), l(5d) 
7(6d), 2(7d) 

151,916 

2.45 

0.069 

.801 * 0.004 


(1) k(ih)=k arcs of £ hours 
k d) = k arcs of i days 

( 2 ) a posteriori variance of unit weight 


For both of these baselines the pass-baseline geometry is weak and 
that of baseline 7110 - 7105 is even worse since passes intersecting this 
baseline at *20 to *60 degrees are only possible, and many of them are 
located close to the perpendicular bisector plane of this baseline. 

For both of these baselines steady state response was not possible 
through a short arc solution, therefore, long arc solutions were 
performed on the basis of 42, 54, 62 and 72 passes and with maximum 
arc lengths up to two, three, and seven days respectively. 

Examination of the results listed in rows 2-4 of Table 16 reveals 
that with maximum arc length up to two days steady state response has 
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been reached when 54 passes are available and the a posteriori standard 
deviation of the estimated baseline is 3 mm. Increasing the maximum arc 
length up to three days leads, as expected, to a faster steady state 
response since such a response is close to being reached with 42 rather 
than with 54 passes (Table 16, rows 2 and 5). When maximum arc length 
of 7 days is allowed, we see that steady state response has clearly been 
reached on the basis of 42 passes (Table 16 cont'd, row 1). However, 
the a posteriori standard deviation is greater than one, thereby 
suggesting that either accumulated orbital errors and/or possible 
ill-conditioning may be on their way up to corrupt the solution. 
Ill-conditioning, however, seems to be at work since many components of 
the recovered initial state vectors exhibit correlation at the .9999 level 
resulting mainly from the weak geometry and the inhomogeneity of the 
arc lengths employed in the solution shown in the 1st row of Table 16 
cont'd. This is also confirmed, when on the basis of 54 passes (2nd 
row), the estimated baseline length jumps by 11 cm and the a 
posteriori standard deviation of unit weight jumps to 3.33. These 
behavior is caused by the high correlations existing among the 
components of the initial state vectors of those short or long arcs, the 
geometry of which became worse by the additional long arcs when 54 
rather than 42 passes were employed (Table 16 cont’d, rows 1-2). 
Inclusion of additional observations, however, seems to lead the response 
of this solution to the right direction, as it is seen from the 2nd 
through 4th row of this table, but at a very slow pace, since the 
geometry is still too weak for a steady state response to take place. 

The divergent response resulting from weak geometry manifests itself 
when baseline 7110 - 7105 is estimated on the basis of 56, 60 and 65 
passes and with maximum arc lengths up to three, five and seven days 
respectively. For many of these passes the geometry with respect to 
baseline 7110 - 7105 is very weak (Section 4.4.1). This in turn leads to 
high correlations among the components of the corresponding initial 
state vectors. With arc lengths up to three days, steady state response 
seems to have been reached because with the incorporation of additional 
observations the length of the estimated baseline just changes at the 
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Table 17 Steady State Response of Baseline 7110 ~ 7105, Long-Arc 
Mode, Passes Within *20* - *60* 

All length units in meters. 


No. of Integration 
Passes Lengths ( * ) 

No. of 
Observ. 

(2) 

/r2 

<r 0 

RMS of All 
Residuals 

SRD Baseline 
3559743. 

Maximum Arc Length ~ 

3 days 




56 12(lh), l(5h) 

l(8h), 5(ld) 
3(2d), 4(3d) 

120,111 

0.58 

0.036 

.594 * 0.003 

60 12(lh), l(5h) 

l(8h), 4(ld) 
3(2d), 5(3d) 

136,356 

0.56 

0.036 

.615 * 0.003 

65 15(lh), l(5h) 

l(8h), 5(ld) 
3(2d), 5(3d) 

146,956 

0.55 

0.035 

.619 ± 0.003 

Maximum Arc Length = 

5 days 




56 ll(lh), l(5h) 
l(8h), 3(ld) 
2(2d), l(3d) 
3(4d), l(5d) 

120,114 

0.64 

0.038 

.542 * 0.003 

60 ll(lh), l(5h) 
l(8h), 2(ld) 
2(2d), 2(3d) 
3(4d), l(5d) 

136,359 

0.63 

0.038 

.554 * 0.003 

Maximum Arc Length = 

7 days 




56 8(lh), l(ld) 

l(2d), 2(4d) 
5(6d), 2(7d) 

120,118 

1.25 

0.053 

.637 * 0.004 

60 8(lh), l(ld) 

2(2d), 2(4d) 
5(6d), 2(7d) 

136,362 

1.14 

0.051 

.615 ± 0.004 


(1) k(ih)=k arcs of i hours 
k (i d) = k arcs of ^ days 

( 2 ) a posteriori variance of unit weight 


centimeter level and the associated standard deviations assume 3 mm 
values (Table 17, rows 1-3). However, correlations at the 0.9999 level 
exist among the components of the initial state vectors of those arcs 
that are primarily composed of passes crossing this baseline close to its 


160 


perpendicular bisector plane. These correlations become larger when 
maximum arc lengths up five and up to seven days are allowed and 
when additional observations result in sharp inhomogeneity in regard to 
the arc lengths incorporated in the same solution (Table 17, rows 4-5, 
and 6-7). Thus, on the basis of 65 passes and when arcs up to five or 
seven days are allowed, the existence of high correlations leads to a 
singular normal equation matrix, thereby implying divergent response. 
Because of this, solutions were not possible with those arc lengths and 
therefore they are not shown in Table 17. 

Thus, it is questionable whether for long baselines (>3500 km) 
steady state response through single SRD baseline solutions can be 
reached. For shorter baselines, however, steady state response can be 
reached even if passes parallel to the baseline do not exist. 

The ill-conditioning affecting longer baselines could be prevented in 
a network solution if there exist observations from baselines that are 
parallel to the passes responsible for the ill-conditioning. Such 
solutions, however, for the reasons mentioned in the previous section 
could not be performed in the course of this study. 

4.5 Baseline Ck>mparison 

The baseline lengths estimated with the geometric method are 
independent of any orbital errors and inconsistencies affecting the 
implementation of the terresrial reference frames; therefore, these 
lengths (Section 4.3.2) will constitute the standards of comparison when 
assessment is made about the accuracy of the baselines estimated via 
the SRD method. Both the geometric and the SRD baseline estimates will 
also be compared with those obtained via long arc range dynamic 
methods, on the basis of the MERIT data set, by the Central Institute 
for Physics of the Earth (ZIPE) in East Germany and by the Center for 
Space Research (CSR) at The University of Texas. 

Both of these centers used in their solutions the gravitational 
constant proposed by MERIT standards (i.e., GM =3.98600448 * 10 

m’s“*) which is different from that employed in the single SRD baseline 
solutions (Section 2. 2. 5.1). A change in the gravitational constant will 
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not only affect the metric scale (the semimajor axis a) but also the 
dynamic scale (mean anomalistic motion n) according to Kepler’s modified 
third law : 

n^a^ =GM(1+X) (4-9) 

where X is a small parameter depending on the satellite orbit. It is 
uncertain how the effect of this change will be divided between the 
metric and the dynamic scales. For regional baseline estimation 
however, the metric scale is primarily controlled by the velocity of light 
which is the same for all the SRD, ZIPE and (CSR)85L01 solutions. For 
this reason, changing the value of the gravitational constant (GM) in the 
SRD solutions presented in the previous section in accordance with 
MERIT Standards will affect at the centimeter level the lengths of only 
the two very long baselines (i.e., 7109-7105 and 7110-7105). The 
remaining baselines are affected either at the few millimeter or 
submillimeter level depending on the length of the estimated baseline 
and on the duration of the arcs employed in the corresponding solution. 
Nevertheless, for the sake of consistency the baseline differences shown 
in Table 18 have been adjusted to the same value of the gravitational 
constant by reestimating the SRD baselines on the basis of the (GM) 
value proposed by the MERIT Standards. 

Table 18 lists the differences of those baselines only for which 
steady state response has been reached either through the SRD or the 
geometric solutions. The first column of this table lists the baseline 
id’s, the second one contains the length of the baseline rounded to the 
nearest meter, the third and fourth columns contain the differences of 
the baselines estimated through the SRD method and through the range 
dynamic method by both ZIPE and GSR, the fifth and sixth columns list 
the differences of the geometric and dynamic ZIPE and GSR baselines. 
The last column contains the differences of the baselines obtained 
through the SRD and the geometric methods. The baseline differences 
not listed in this table were not computed because steady state response 
for the SRD and/or the geometric method was not possible or because 
one of the stations constituting a baseline endpoint was not included in 
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Table 18 Baseline Differences 


Baseline 

Length (m) 
to nearest 
meter 

ZIPE 

-SRD 

(CSR)85L01 

-SRD 

ZIPE 

-GEOM 

(CSR)85L01 

-GEOM 

SRD 

-GEOM 

7109-7110 

883602. 

-0.03 

-0.03 

-0.03 

-0.03 

0.00 

7109-7265 

627044. 

— 

— 

0.00 

0.04 

— 

7109-7886 

8. 

-0.01 

0.00 

-0.02 

-0.01 

-0.01 

7110-7122 

1437139. 

-0.03 

-0.03 

-0.02 

1 

O 

o 

to 

0.01 

7110-7220 

15. 

— 

-0.04 

— 

-0.01 

0.03 

7110-7265 

274069. 

-0.04 

-0.07 

-0.03 

-0.06 

0,01 

7122-7265 

1663983. 

— 

— 

-0.02 

-0.06 

— 

7122-7886 

2280720. 

— 

— 

0.01 

0.02 

— 

7220-7265 

274066. 

— 

— 

— 

0.00 

— 

7265-7886 

627049. 

— 

— 

0.01 

0.06 

— 

7110-7886 

883606. 

-0.01 

0.00 

-0.01 

0.00 

0.00 

7110-7086 

1198291. 

0.01 

0.02 

— 

— 

— 

7109-7105 

3703352. 

0.06 

0.05 

— 

— 

— 


the corresponding dynamic method. For instance, the difference (ZIPE - 
3RD) for baseline 7220-7265 could not be computed because steady state 
response for this baseline was not possible through a single 3RD 
baseline solution; also for the same baseline the difference (ZIPE - 
GEOM) is not included in Table 18 because station 7220 was not listed in 
the reported ZIPE solutions (Montag et al., 1985), 

The baseline differences between the 3RD, ZIPE and GSR solutions 
shown in the third and fourth columns are negative for north-south 
baselines (rows 1-11) and positive for the east-west ones (rows 12-13) 
(Fig. 18). 3ince some of these differences are larger in magnitude for 
shorter baselines, scale difference between 3RD and CSR or ZIPE 
solutions would account for part of these differences. The remaining 
differences, at the 2 cm or 3 cm level, must be caused by orbital errors 
affecting mostly the dynamic solutions, because the 3RD solutions 
(column 7) are clearly closer to the geometric solution than both ZIPE 
and C3R solutions (Table 18, columns 5-6). The large differences of 6 
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cm between (CSR)85L01 and geometric solutions (column 6) are associated 
with station 7265 (MOHAVE). This station during the MERIT Main 
Campaign was equipped with a TLRS-1 laser instrument which 
experienced many problems, as was confirmed in the present study when 
editing the data with the data snooping procedure (Section 3.5). 
Therefore, it is very likely that erroneous observations from this station 
are still present in the (CSR)85L01 solution. This in turn would explain 
these relatively large baseline differences. 

Although the geometric solution (Table 8, Type C3) used in Table 18 
has been overconstrained to the (CSR)85L01 solution, SRD baseline 
estimates are on the average the "best" ones as compared to those of 
either the ZIPE or the (CSR)85L0l solutions. The term "best" indicates 
that a solution is the one closest to a geometric solution at its steady 
state. However, ZIPE baseline estimates when compared to those 
obtained in the (CSR)85L01 solution, are closer to the geometric solutions 
and therefore more accurate. 

The root mean square of the differences between the geometric 
baseline estimates and those of the SRD, ZIPE and CSR are 1 cm, 2 cm, 
and 4 cm respectively. 

Since these differences are based on baselines up to 1500 km, it is 
fair to state that SRD baseline estimates are at least as accurate as 
those obtained through the range dynamic mode methods. 

4.6 RESPONSE OF THE SRD METHOD TO THE SIMPLIFICATIONS OF THE 

ORBITAL MODEL 

It was mentioned in Section 2.2.5 that the aim of the present study is 
not to estimate Lageos’ orbit with the highest degree of accuracy but 
rather to employ models as simple as possible and yet be able to 
recover baselines with an accuracy compatible to that of the 
observations. 

Since the temporal variations of the baseline endpoints have been 
accounted for to the required degree of accuracy (Section 3.3), and 
since any inconsistencies in the implementation of the Terrestrial 
Reference Frame do not affect the SRD observables (Section 4.4.2), the 
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errors affecting the baselines estimated via the steady state response of 
the SRD method (Section 4.4.2) result only from those orbital errors 
accumulated over the integration periods and not cancelled out in the 
computation of the SRD observable. Thus, the question that should be 
addressed and investigated is twofold: 

1) Is the sophistication of the orbital model, employed in the present 
study (Section 2.2.5) sufficient to result in baselines the accuracy of 
which is compatible to that of the observations? 

2) If the answer is yes, how much could the sophistication of the orbital 
model be reduced without affecting the accuracy of the estimated 
baselines? If the answer is no, how much should the sophistication 
of the orbital model be enhanced so the estimated baselines have an 
accuracy compatible to that of the observations? 

The errors affecting the estimated baselines are propagated from 
those affecting the corresponding Cartesian coordinates. These errors, 
which originated by the erroneous constraints imposed on a large 
number of estimable quantities entering in the computation of the SRD 
observable (eq. 2-60), tend to accumulate as the employed integration 
periods become longer. 

For baselines of moderate length (<2000 km) accumulated radial 
errors are cancelled out almost totally in the computation of the SRD 
observable. These errors are propagated almost unaltered into the 
computed range observable (Pavlis, 1982). Depending, however, on the 
location of the observed satellite positions, accumulated latitudinal and 
longitudinal errors may affect the computed value of the SRD observable 
worse than they affect the range observable (ibid.) For shorter 
baselines (<200 km) the computed SRD observable is less affected by all 
these three errors (i.e., radial, latitudinal and longitudinal). 
Consequently, the answers to the questions posed in the beginning of 
this section depend on the length of the estimated baseline, on the 
magnitude of the accumulated orbital errors, and on whether these 
errors are radial, latitudinal or longitudinal. If these errors are 
primarily latitudinal and longitudinal, these answers depend also on the 
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relative orientation of the estimated baselines and its observed passes. 
This is the result of the anisotropy characterizing the latitudinal and 
longitudinal error surfaces. The orientation of the estimated baseline is 
not important for the propagation of the radial orbital errors because 
the radial error surfaces are isotropic. Since all of the baselines of 
Table 19 are extended in both north-south and east- west directions (Fig. 
18), they will be affected by all three orbital errors (radial, latitudinal 
and longitudinal). 

The estimable quantities, entering the satellite perturbations and 
being neglected in the simplification process, implicitly assume zero 
values in the resulting orbital model. If these constraints result in 
radial errors, then depending on the length of the estimated baselines, 
these errors should be relatively large in order to affect those baselines 
beyond the centimeter level. This in turn would allow longer lengths of 
continuous integration and, therefore, a reduced number of observations 
would be necessary to achieve steady state response of the SRD method 
(Section 4.4.2). However, if these constraints result, in addition to 
radial errors in both latitudinal and longitudinal errors, then even 
smaller lengths of continuous integration may affect baselines of 
moderate length (<2000 km) beyond the centimeter level. To set up the 
guidelines as to what simplifications of the orbital model can be applied 
without affecting the accuracy of the estimated baseline beyond the 
centimeter level, several tests have been performed the results of which 
are shown in Table 19. 

Table 19 contains the baseline differences obtained as the orbital 
model was simplified from one containing a 12x12 gravity field, the 
direct point mass (PM) effects of Sun and Moon, the tidal(TD) effects 
due to Sun and Moon, the solar radiation (SR) pressure effects and the 
along-track (AT) acceleration effects, to that containing only a 2x2 
gravity field and the direct (PM) effects of the Sun and the Moon. The 
first column of Table 19 lists the orbital models employed to estimate the 
baselines which subsequently were differenced from those estimated on 
the basis of orbital model shown on the title of this table (i.e., 
12xl2+(D). 
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Table 19 Baseline Differences (in meters) With Respect to Those 
Computed Using an Orbital Model Including a 12x12 
Gravity Field + (1) 


Force Model Gravity 

7110-7265 

7109-7110 

7110-7122 

Gravity Field + ( ) 

(274069.50) 

(883602.25) 

(1437139.30) 

12x12 + 

(2) 

0.01 

0.00 

0.00 

12x12 + 

(3) 

0.02 

0.00 

0.04 

12x12 + 

(4) 

-0.08 

0.00 

0.10 

12x12 


5.14 

-0.02 

-3.15 

10x10 + 

(2) 

0.01 

0.00 

0.01 

10x10 + 

(4) 

— 

0.00 

— 

8x8 + 

(2) 

0.03 

0.00 

0.01 

8x8 + 

(4) 

— 

0.00 

— 

6x6 + 

(2) 

1.02 

0.00 

2.02 

6x6 + 

(4) 

— 

0.01 

— 

4x4 + 

(2) 

0.90 

0.00 

3.37 

4x4 + 

(4) 

— 

-0.01 

— 

3x3 + 

(4) 

— 

-0.02 

— 

2x2 + 

(4) 

— 

-0.07 

— 


(1) 

(PM) 

+ 

(TD) 

+ (SR) + (AT) 

PM 

(2) 

(PM) 

+ 

(TD) 

+ (SR) 

TD 

(3) 

(PM) 

+ 

(TD) 


SR 

(4) 

(PM) 




AT 


= point mass effects of sun & moon 
= tidal effects due to sun & moon 
= solar radiation pressure effects 
= along-track acceleration effects 


The resulting differences are shown in the corresponding rows for only 
three baselines estimated on the basis of arcs up to seven days (Table 
11, row 4), up to one hour (Table 12, row 4) and up to three days 
(Table 15, row 4). 

Elimination of the AT acceleration, the SR pressure, the TD 
acceleration and the direct PM effects of the Sun and the Moon, 
introduces into the resulting orbital errors secular, long-period and 
short-period terms. Short-period terms directly related to Lageos' mean 
anomaly are common in all of these errors. 

Since the eccentricity of Lageos' orbit is very small ("0.0039), 
elimination of the AT acceleration results primarily in latitudinal and 
longitudinal orbital errors (Obenson, 1970). The long periods of the 
resulting errors range from 66 days to 1100 days (Smith et al., 1985). 
Elimination of the SR pressure acceleration results also in radial, 
latitudinal and longitudinal orbital errors. Besides the additional 
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short-period terms related to the "ramp-like" behavior of the solar 
radiation pressure, the long periods of the resulting orbital errors are 
directly related to the motion of the Earth around the Sun and to the 
motion of the ecliptic in space (Musen, 1960; Kozai, 1961; El’Yasberg, 
1967; Blitzer, 1970)* In addition to the short-period perturbations 
related to the daily motion of the Earth, the nonstationary disturbances 
of the Earth's potential due to the attraction of the Sun and the Moon 
perturb the satellite periods with periods greater than a week (Kozai, 
1973; Goad, 1977)* As a result radial, latitudinal and longitudinal orbital 
errors, having those periods will be introduced from the elimination of 
the TD acceleration. Finally, elimination of the direct PM effects of the 
Sun and the Moon introduces again the three types of orbital errors 
having, besides the short periods associated with Lageos' mean anomaly, 
intermediate periods associated with the motion of the Earth around the 
Sun(^ multiple of 180 days) and the motion of the Moon around the 
Earth('^ multiple of 14 days), and long periods associated with the space 
motions of the ecliptic and the orbital plane of the Moon (Kozai, 1959; 
El'Yasberg, 1967; Fisher, 1971; Blitzer, 1970). 

Therefore, elimination of each of the mentioned accelerations 
introduces radial, latitudinal and longitudinal orbital errors. Since 
periods longer than a week are present, these errors will accumulate as 
the integration periods increase from one hour to seven days. The 
effects of these errors on the estimated baselines are shown in Table 19 
(rows 1-4). 

The perturbations caused in the motion of the satellites by the 
ocean tides can reach as much as 20% of the perturbations caused by 
the tides of solid Earth (Musen, 1973). Inspection of Table 19, row 3, 
reveals that elimination of ocean tidal effects from the orbital model, as 
is the case in this study (Section 2.2.5), will hardly affect the estimated 
baselines at the centimeter level. 

Therefore, centimeter level accuracy for baselines up to 1500 km, 
estimated via long-arc solutions, can be achieved if in addition to 
gravitational effects of the Earth, the SR pressure, the TD and the 
direct PM effects of Sun and Moon are included in the orbital model. 
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However, if steady state response is possible via a short arc solution, 
centimeter level accuracy for baselines up to 1000 km can be achieved if 
in addition to the gravitational effects of the Earth only the direct PM 
effects of the Sun and the Moon are included (Table 19, column 3). 

The question as to how much the gravity field of the Earth can be 
reduced without affecting at the centimeter level the estimated baselines 
is investigated by performing several solutions, the results of which are 
shown Table 19 (rows 5-14). In these solutions the gravity field is 
being reduced by two degrees at a time while the SR pressure, the TD 
and the direct PM effects of the Sun and the Moon are included in the 
orbital model (Table 19, rows 5, 7, 9 and 11). For baseline 7109-7110, 
estimated via short-arc solutions, an additional solution has been 
performed with each reduced gravity field and with only the direct PM 
effects of the Sun and the Moon (Table 19, column 3, rows 6, 8, 10, 12, 
13 and 14). Eliminating two degrees at the time introduces secular 
orbital errors due to elimination of the even zonal harmonics, 
long-period errors due to the elimination of the zonals, and short 
periods due to the elimination of all of the harmonics included in those 
degrees. Thus, as the length of continuous integration increases radial, 
latitudinal and longitudinal orbital errors tend to increase to the extent 
that the estimated baselines will be affected beyond the centimeter level. 

Careful study of the results shown in Table 19 reveals that 
baselines of up to 1500 km estimated via the SRD method will be affected 
at just the centimeter level if the orbital model includes : 

1) in short arc solutions: 

gravity field 4x4, and 

the direct point mass effects of the Sun and the Moon. 

2) in long arc solutions with arcs up to three days: 

gravity field (8x8), 

the PM effects of the Sun and Moon, 

TD effects, and 
SR pressure effects. 
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3) in long-arc solutions with arcs up to seven days: 
gravity field (10x10), 
the PM effects of the Sun and Moon, 

TD effects, and 
SR pressure effects. 

Therefore, the sophistication of the orbital model employed in the 
present study results in baselines having centimeter-level accuracy. 
This accuracy is compatible with the accuracy of the laser observations 
employed in this study (Section, 4.2). This constitutes the answer to 
question (1); question (2) has already been answered by the conclusions 
stated above. 
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Chapter 5 


CONCLUSIONS AND RECOMMENDATIONS 


5.1 CONCLUSIONS 

The severe requirements of the geometric solutions to have simultaneous 
observations from four or more stations (Section 2.1.3), results in 
rejection of a large number of nonsimultaneous observations. The 
rejection is even greater with weather dependent satellite observations 
not specifically designed for simultaneous tracking as happens with the 
laser range observations. 

Although during the MERIT Main Campaign many stations committed 
themselves to collecting simultaneous observations (Section 3.4), these 
observations were not enough to yield strong geometric solutions 
because of their sensitivity to configurations being close to those 
leading to singularities B and C (Sections 2.1.3 and 4.3.1). Because of 
these configurations, a steady state response in minimum constraint 
geometric solutions was possible only for two very short baselines. 
However, stronger implementation of the scale through an overconstraint 
solution (Table 8, solution C3) resulted in a steady state response for 
baselines of up to 2280 km (Section 4.3.2). These baselines formed the 
standards of comparison in assessing the accuracy of those obtained via 
both the SRD semidynamic and range dynamic mode methods (Section 
4.5). Although the scale in this overconstraint solution was partly 
implied by that of the (CSR)85L01 (Section 4.3.2), it is feasible in the 
forseable future to incorporate the VLBI scale on the basis of only one 
or two baselines. This in turn could lead through a geometric solution 
to the estimation of a large number of baselines of compatible accuracy. 
These baselines could be effectively used to assess the accuracy of 
those obtained via the range dynamic methods. This practice might be of 
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great importance with the upcoming millimeter accuracy laser range 
systems. 

As a result, regular baseline estimation through geometric soliitions 
using laser range observations is impossible and therefore, since 
geometric solutions are free of errors affecting the orbit and the TRF 
frame they should be used, whenever possible, as standards of 
comparison. 

In contrast to the geometric method, the SRD semidynamic mode 
method can be effectively used for regional baseline estimation, 
especially with laser range observations to Lageos, since the altitude of 
its orbit makes it possible for stations as far as 3703 km apart to collect 
enough simultaneous observations to result in a steady state response of 
the SRD method. The number of observations required for such a 
response to take place is a function of 

the baseline length, 

the geometry and the length of the arcs employed, 
and the accuracy of the available observations. 

Table 20 shows the number of passes used in the present study to 
achieve steady state response for different baseline lengths. In 
general, as the length of the estimated baseline increases, more arcs, 
stronger geometry, and more observations are required for a steady 
state response to take place. The speed of the steady state response 
will primarily depend on the number, distribution and accuracy of the 
available observations, and on the length of the arcs employed in the 
SRD solutions. 

For instance, when the geometry is favorable (i.e., parallel passes 
are available) and the observations of the baseline end stations have an 
accuracy of about 3 cm, baselines of 883 km can be estimated via both 
short arc and long arc solutions (Table 12). This is possible on the 
basis of only 12 passes collected within a period of a week. However, 
as the accuracy of the observations decreases, even in the light of 
favorable geometry, the steady state response of a baseline having 
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Table 20 Steady State Response of the SRD Method 


No. of Passes 
Needed 

Approx. Occupation Time 
in This Dataset 

steady State 
Up to km 

10 - 15 

1 week 

1000 

20 - 25 

3 months 

1500 

25 - 30 

3 months 

2500 

50 - 55 

8 months 

3500 


about the same length as that in Table 12 (see Table 13) was not 
possible on the basis of approximately the same number of observations 
and through a short arc solution because the accuracy of its end 
stations was worse than that of Table 12 (see Table 2), Nevertheless, a 
steady state response for this baseline was possible through a solution 
employing arcs of up to two days long. For longer baselines of up to 
1437 km, the end stations of which recorded observations of 3 cm and 9 
cm accuracy respectively, a steady state response in the light of 
favorable geometry was possible on the basis of 22 passes and only 
through a long arc solution employing arcs up to three days long (Table 
15). 

when the baseline-pass geometry is not favorable, (i.e., passes *30* 
to *50*) and the accuracy of the observations of the baseline end 
stations is approximately 3 cm and 8 cm respectively, a steady state 
response for a 274 km baseline was possible on the basis of 17 passes 
and only through a long arc solution employing arcs up to four days 
long (Table 11). For another baseline of about 1198 km having 
unfavorable pass-baseline geometry, a steady state response was 
possible on the basis of 40 passes and through a solution employing 
arcs up to seven days long (Table 14). The number of observations 
required for the steady state response of this baseline is substantially 
reduced as compared to those required for all the other baselines shown 
in Tables 9 through 17. This is not surprising because the stronger 
geometry implied by the 40 passes compensated for the reduced number 
of observations. However, if more observations could have been 
collected during each pass, fewer passes would be required for a steady 


173 




state response. 

For very long baselines of about 3703 km, although passes of 
favorable geometry are not available, a steady state response was still 
possible on the basis of 54 passes and through a solution employing 
arcs up to one day long (Table 16). This response, however, started 
diverging when for the determination of this baseline arcs up to seven 
days long were employed (Table 16). Furthermore, for a baseline of 
about 3559 km, the pass-baseline geometry of which is very weak (i.e., 
passes within *20* to *60*), a steady state response was possible only 
with arcs up to three days (Table 17). Therefore, for so long baselines 
and when the geometry is fairly weak care should be exercised as to 
whether or not a steady state response is possible. 

From the above discussion it is evident that for baselines up to 
1500 km, steady state response can be reached without any problems 
even if the geometry is not favorable. However, tracking passes of 

favorable geometry with accurate (i.e., 3 cm) and high repeatability laser 
instruments increases the resolution of baseline recovery to weekly 

estimates. Although this was shown only for baselines of up to 883 km, 

baselines of up to 1437 km can also be recovered on a weekly basis, 

because steady state response for these baselines was reached on the 
basis of about 69,000 observations (Table 15). 

In addition to the increased temporal resolution the comparison of 
the 3RD method with the geometric solution (Table 8, solution C3) 
evidently demonstrated an accuracy of 1cm (Table 18), accomplished on 
the basis of a simple modeling and a limited orbit adjustment (Section 
4.6). This demonstrates the potential of the 3RD method for accurate 
differential positioning, and the insensitivity of the 3RD observable not 
only to the errors affecting the orbit but also to those affecting the 
implementation of the TRF frame. The simple modeling of the orbit and 
its limited adjustment result in relaxing requirements in regard to 
availability of sophisticated software and extensive computer facilities 
which in turn makes it feasible to implement this method in the PC 
environment. 
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The stability characteristics of the SRD single baseline semidynamic 
solutions make it possible to estimate baselines of up to 1437 km without 
any need for a local support tracking network which is required in 
regional baseline estimation via the range semidynamic mode method 
(Christodoulidis et al., 1981). 

The steady state response of baseline 7109-7105 (Table 16), reached 
on the basis of weak pass-baseline geometry strongly suggests that with 
parallel passes single SRD baseline solutions may very well lead to 
accuracy at the 2 to 3 cm level for baselines of up to 3500 km or 4000 
km. 

The stability characteristics and the high accuracy of the SRD 
method make it ideal for regional baseline estimation. This in turn 
makes it possible to obtain baselines up to 1500 km on a weekly basis 
free of the fluctuations usually affecting the monthly estimates of the 
range dynamic method. The stability characteristics could be greatly 
enhanced with a network solution which on the basis of SRD normal 
points would also substantially reduce the bulk of the compiitations 
without affecting the speed of the steady state response since the 
almost noise free SRD normal points will lead to a faster steady state 
(Section 4.4). However, with a network solution there might exist, 
baselines for which simultaneous observations will not be available or 
the SRD observables will contain very little information for those 
baselines. This happens when the intersection of the visibility regions 
of the baseline end points is either an empty set or one having very 
few points. If a steady state response for those baselines has been 
reached it will entirely depend on the strength of the network which in 
the semidynamic environment tends to depend on the geometric strength 
of the observations which may not lead to a steady state response for 
longer baselines as it was the case with the geometric solutions. Steady 
state response of those baselines would require increasing the 
integration period. If this increase leads to a steady state, the 

estimated lengths for those baselines would primarily depend on the 
orbital strength, and therefore they would be sxjsceptible not only to 
accumulated orbital errors but also to any inconsistencies affecting the 
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implementation of the TRF frame. However, employing the same orbital 
model as in the range dynamic solutions, the SRD baseline estimates will 
be less affected by the accumulated orbital errors. 

In a network solution care should be exercised to account for the 
correlations resulting from the redundancy of having pairs of stations 
with common observations. 

In the implementation of the SRD method a great deal of effort was 
spent to edit the data and create the SRD observables. Editing the 
observed laser ranges before generating those observables is crucial 
because generating them on the basis of erroneous laser ranges leads to 
the distribution of these errors among all of the generated SRD 
observables and therefore it will be difficult, if not impossible, to 
effectively edit the SRD observables in the final orbit adjustment. 

The effects of gaps should also be well controlled and should be 
kept below the noise level of the observations. This was the reason 
why so much care was taken in the generation of SRD observables 
(Chapter 3). This in turn increased considerably the bulk of the 
computations required for the implementation of the SRD method. 

These two factors do not by any means limit the potential of the 
SRD method because through proper arrangements either full rate SRD 
observables or preferably SRD normal points could be very easily 
generated with some slight modifications of the software employed by 
the computational centers responsible for the generation of the range 
normal points. With SRD normal points at our disposal a network SRD 
solution could be performed even in a PC environment, and since normal 
points are almost noise free they would lead to a faster steady state 
response of the SRD method. 


5.2 RECOMMENDATIONS 

Since the SRD method, based on a relatively simple orbital modeling and 
a limited orbit adjustment, is for regional baseline estimations at least as 
accurate as the dynamic mode methods, it is recommended that for such 
applications serious consideration should be given to regularly implement 
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this method for accurate differential positioning. If such a step is to 
be taken it is also suggested that: 

Further effort should be made by the observing stations to 
achieve simultaneous tracking. 

The direction of the baselines being estimated should be chosen, if 
this is feasible, to closely resemble the two main Lageos groundtrack 
directions. 

SRD normal points should be generated by the computational center 
responsible for the generation of range normal points. 

Research the response of the method, on the basis of normal points, 
when different network configurations are employed. 

Projects specifically designed for studies of regional plate tectonic 
motions, determined on the basis of laser range observations to Lageos, 
are ideal for implementation of the SRD method. Such a project is 
currently under way in the area around the Mediterranean Sea 
(MEDLAS/WEGENER project), and therefore for this project the SRD 
method could offer an accurate and inexpensive alternative to study 
regional plate tectonic motions. 
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Appendix A 


PARTIAL DERIVATIVES OP THE SRD OBSERVABLE 


Expressing the SRD observable relative to an earth-fixed frame takes 
the following form: 

6pj = (EjJ-Eja)^ - (EjiT-Ep)’i (A-1) 

with 

Ejj = Xj - Ya and Eji = Xj - Yi (A-2) 

where Xj, Ya and Yt are the earth-fixed position vectors of the satellite 
at the epoch j and stations (2) and (1) respectively. 

Differentiating 5/>j with respect to Ya one obtains: 

d6pj 9<5pj dEja 1 _ Eja"^ 

SY2 ~ Eja aY2 ~ ^ ' Pj2 

Thus 

36p. (Xj-Ya)T 

= (A-3) 

3Ya Pj2 

where pja is the distance from station (2) to satellite position J. 
Suppose that 

Ro = (Xo, Yo, Zo) (A-4) 

and 

Ro = (Xo, Yo, Zo)^ (A-5) 

are the initial position and velocity vectors with respect to an inertial 

frame* Differentiating equation (A-1) with respect to the initial state 
^ • 

vector (Ho, Ro) one obtains: 
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(A-6) 


3<5pj 

3(Ro,Ro) 3 Ej 2 aXj 3(Ro,Ro) ^ ^Ej/jXj ’T^l^X) 


From equ. (A-1) and (A-2) it follows that: 

35pj 1 _ XjT-YaT 

3EJ2 PJ2 Pj2 

3Ej 2 

3^ ' ^ 


3Xj 

•3Xj 

3Xj] 

(3(SNP)Rj 3(SNP)Rj] 

3(Ho,Ro) ' 

3Xo’ 

Oj 

N* 

o 

1 

1 aXo ’ azo J 


(A-7) 


(A-8) 


(A-9) 


where pja is the distance from the station (2) to the satellite position j; 
Rj is the inertial position vector of the satellite at the epoch j; and S, 
N» P are the earth rotationi nutation and precession matrices 
respectively. Similarly to eq. (A-7), the partial 36pj/3Ej, takes the 
following form: 


36 pj XjT-Y.t 
3B j 1 P j 1 


(A-10) 


where pj, is the distance from station (1) to satellite position j. 

Substituting eq. (A-7) - (A-10) into eq. (A-6) the following formula 
results: 


35pj 

f(Xj-Y2)T 

(Xj-’'-)'] 


3Xj 3Xj 

3(Ro,Ro) ' 1 

[ Pj2 



aXo’ azo 


(A-11) 


or 

3(Ro,Ro) 


(Xj-Y2)T 

Pj2 


(Xj-Yx)T 


Pji 


(SNP) 


3Xo’ 


3Rj 

. (SNP)^ 
BZo 


(A-12) 


Thus, from eq. (A-11) and (A-12) it is easily deduced that the partial 
derivative of the SRD observable with respect to the first coordinate 
(i.e., Xo) of the initial state vector (Ro.Ro) takes the following form: 
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(A-13) 


3«Pj 

[Xj-Y, 

Xj-Yx] 

T 3Rj 

, /cvrp\ 

aXo " 

Pja 

Pji J 

^^^^3Xo 

36pj 

[Xj-Y, 

Xj-Yxl 

T 3Xj 

axT* “ 

Pj2 

Pji J 

3Xo 


Similarly one can obtain the partial derivative of the SRD observable 
with respect to all of the remaining componenets of the inertial initial 
state vector. 
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Appendix B 


SENSITIVITY OP THE PARTIAL DERIVATIVES WITH RESPECT TO THE 

INITIAL STATE VECTORS 


With the tendency of the observed ranges pja and pji from stations (1) 
and (2) to become equal (Section 4.4.1). 


PJ2 “ Pjl - Pj 
equation (A-11) becomes 

36pj (Yi-Ya)T pXj 

3(Ro,Ro) ' Pj aZoJ 


(B-1) 


(B-2) 


This equation can also be expressed as follows 


^ (®o » ®o ) 


1 

— (Y»-Ya)-(Di{o. .... H'zi 


where (•) and (Dx^, .... D^o) denote the dot product and the partial 

derivative vectors • • • • at the epoch j respectively. 

'■•'Ao aZo-' 

If «xj’» •••» designate the angles between the baseline vector 

and the partial derivative vectors respectively, equation (B-3) takes the 
following form 


afipj 1 . . _ , 

— = — ^ I Yi-Ya I ( IDx^l -cos(wx^) , IDj^lcoswzi) (B-4) 

3 (Ro,Ro) Pi 


Therefore, from epoch to epoch the variations among the entries of 
the columns, corresponding to the components of an initial state vector, 
are controlled by the variations of the projection of the partial deriv- 
ative vectors ®x„, ..., in the direction of the estimated baseline. 
These variations tend to be reduced as the observed part of the pass 


c ^ 3 
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becomes short and/or perpendicular to the estimated baseline. This 
geometry is closer to reality for longer baselines because as the 
baseline length increases the intersection of the visibility regions of the 
baseline endpoints tends to become smaller and closer to the 
perpendicular bisector plane of the baseline, thereby leading to a 
geometry close to also fulfilling assumption (B-1). 
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